Entrega grupal II

Aplicando regresión univariante, modelo saturado, regresión con selección de modelos y KNN y árbol en modo regresión a house_prices.csv

Iván González, Qiqi Zhou y Felipe Reyes (Universidad Complutense de Madrid)https://ucm.es
2023-01-15

Objetivo

El objetivo de esta práctica es predecir la variable continua SalePrice a través de distintos métodos de regresión y de algoritmos en modo regresión.

Paquetes necesarios

Necesitaremos los siguientes paquetes:

# Borramos
rm(list = ls())

# Paquetes
library(skimr) # Resumen numérico
library(tidymodels) # Depuración datos
library(tidyverse) # Modelos
library(outliers) # Outliers
library(themis) # Sobremuestreo
library(parallel) # Fase de validación
library(doParallel) # Fase de validación
library(rpart.plot) # Visualización de árboles
library(performance)
library(ggthemes)
library(glue)
library(vip)
library(ggrepel)
library(Rmisc)

Datos

Los datos que usaremos provienen de los dataset house_prices_train.csv y house_prices_test.csv. Para el análisis exploratorio, bindearemos ambos dataset para poder analizar el total de las observaciones.

# Cargamos ambas particiones
house_train <- read_csv(file = "/Users/leztin/Desktop/DATOS/house_prices_train.csv")
house_test <- read_csv(file = "/Users/leztin/Desktop/DATOS/house_prices_test.csv")

# Creamos en test una columna con la variable objetivo `SalePrice`
house_test$SalePrice <- NA

# Bindeamos las particiones en un archivo completo
house_complete <- rbind(house_train, house_test)

Análisis exploratorio inicial (numérico)

Antes de tomar cualquier decisión con los datos, lo primero que haremos será echar un vistazo numérico a cómo se comportan las variables. Dado que dos de los métodos que vamos a utilizar para predecir son el árbol y el KNN en modo regresión, comprobaremos cómo se relacionan nuestras predictoras cuantitativas y cualitativas de cara a la creación de la receta y con tal de poder recategorizarlas.

Variables

Lo primero es conocer nuestras variables.

glimpse(house_complete)
Rows: 2,919
Columns: 52
$ Id            <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,…
$ MSSubClass    <dbl> 60, 20, 60, 70, 60, 50, 20, 60, 50, 190, 20, 6…
$ MSZoning      <chr> "RL", "RL", "RL", "RL", "RL", "RL", "RL", "RL"…
$ LotFrontage   <dbl> 65, 80, 68, 60, 84, 85, 75, NA, 51, 50, 70, 85…
$ LotArea       <dbl> 8450, 9600, 11250, 9550, 14260, 14115, 10084, …
$ Street        <chr> "Pave", "Pave", "Pave", "Pave", "Pave", "Pave"…
$ Alley         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LandContour   <chr> "Lvl", "Lvl", "Lvl", "Lvl", "Lvl", "Lvl", "Lvl…
$ Utilities     <chr> "AllPub", "AllPub", "AllPub", "AllPub", "AllPu…
$ LandSlope     <chr> "Gtl", "Gtl", "Gtl", "Gtl", "Gtl", "Gtl", "Gtl…
$ Neighborhood  <chr> "CollgCr", "Veenker", "CollgCr", "Crawfor", "N…
$ BldgType      <chr> "1Fam", "1Fam", "1Fam", "1Fam", "1Fam", "1Fam"…
$ HouseStyle    <chr> "2Story", "1Story", "2Story", "2Story", "2Stor…
$ OverallCond   <dbl> 5, 8, 5, 5, 5, 5, 5, 6, 5, 6, 5, 5, 6, 5, 5, 8…
$ YearBuilt     <dbl> 2003, 1976, 2001, 1915, 2000, 1993, 2004, 1973…
$ YearRemodAdd  <dbl> 2003, 1976, 2002, 1970, 2000, 1995, 2005, 1973…
$ RoofStyle     <chr> "Gable", "Gable", "Gable", "Gable", "Gable", "…
$ RoofMatl      <chr> "CompShg", "CompShg", "CompShg", "CompShg", "C…
$ Exterior1st   <chr> "VinylSd", "MetalSd", "VinylSd", "Wd Sdng", "V…
$ ExterCond     <chr> "TA", "TA", "TA", "TA", "TA", "TA", "TA", "TA"…
$ BsmtQual      <chr> "Gd", "Gd", "Gd", "TA", "Gd", "Gd", "Ex", "Gd"…
$ BsmtCond      <chr> "TA", "TA", "TA", "Gd", "TA", "TA", "TA", "TA"…
$ TotalBsmtSF   <dbl> 856, 1262, 920, 756, 1145, 796, 1686, 1107, 95…
$ Heating       <chr> "GasA", "GasA", "GasA", "GasA", "GasA", "GasA"…
$ HeatingQC     <chr> "Ex", "Ex", "Ex", "Gd", "Ex", "Ex", "Ex", "Ex"…
$ CentralAir    <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "…
$ Electrical    <chr> "SBrkr", "SBrkr", "SBrkr", "SBrkr", "SBrkr", "…
$ `1stFlrSF`    <dbl> 856, 1262, 920, 961, 1145, 796, 1694, 1107, 10…
$ `2ndFlrSF`    <dbl> 854, 0, 866, 756, 1053, 566, 0, 983, 752, 0, 0…
$ GrLivArea     <dbl> 1710, 1262, 1786, 1717, 2198, 1362, 1694, 2090…
$ FullBath      <dbl> 2, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 3, 1, 2, 1, 1…
$ HalfBath      <dbl> 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0…
$ BedroomAbvGr  <dbl> 3, 3, 3, 3, 4, 1, 3, 3, 2, 2, 3, 4, 2, 3, 2, 2…
$ KitchenQual   <chr> "Gd", "TA", "Gd", "Gd", "Gd", "TA", "Gd", "TA"…
$ TotRmsAbvGrd  <dbl> 8, 6, 6, 7, 9, 5, 7, 7, 8, 5, 5, 11, 4, 7, 5, …
$ Functional    <chr> "Typ", "Typ", "Typ", "Typ", "Typ", "Typ", "Typ…
$ Fireplaces    <dbl> 0, 1, 1, 1, 1, 0, 1, 2, 2, 2, 0, 2, 0, 1, 1, 0…
$ FireplaceQu   <chr> NA, "TA", "TA", "Gd", "TA", NA, "Gd", "TA", "T…
$ GarageYrBlt   <dbl> 2003, 1976, 2001, 1998, 2000, 1993, 2004, 1973…
$ GarageCars    <dbl> 2, 2, 2, 3, 3, 2, 2, 2, 2, 1, 1, 3, 1, 3, 1, 2…
$ GarageArea    <dbl> 548, 460, 608, 642, 836, 480, 636, 484, 468, 2…
$ GarageQual    <chr> "TA", "TA", "TA", "TA", "TA", "TA", "TA", "TA"…
$ PavedDrive    <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "…
$ OpenPorchSF   <dbl> 61, 0, 42, 35, 84, 30, 57, 204, 0, 4, 0, 21, 0…
$ PoolArea      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ PoolQC        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ Fence         <chr> NA, NA, NA, NA, NA, "MnPrv", NA, NA, NA, NA, N…
$ MoSold        <dbl> 2, 5, 9, 2, 12, 10, 8, 11, 4, 1, 2, 7, 9, 8, 5…
$ YrSold        <dbl> 2008, 2007, 2008, 2006, 2008, 2009, 2007, 2009…
$ SaleType      <chr> "WD", "WD", "WD", "WD", "WD", "WD", "WD", "WD"…
$ SaleCondition <chr> "Normal", "Normal", "Normal", "Abnorml", "Norm…
$ SalePrice     <dbl> 208500, 181500, 223500, 140000, 250000, 143000…
Variable Significado Variable Significado
SalePrice precio de venta de la propiedad en dólares (variable objetivo) Heating tipo de calefacción
MSSubClass clase de construcción HeatingQC calidad y condición de la calefacción
MSZoning clasificación general de zonificación CentralAir aire acondicionado central
LotFrontage pies lineales de calle conectados a la propiedad Electrical sistema eléctrico
LotArea tamaño de la propiedad en pies cuadrados 1stFlrSF pies cuadrados del primer piso
Street tipo de acceso a la carretera 2ndFlrSF pies cuadrados del segundo piso
Alley tipo de acceso al callejón LowQualFinSF pies cuadrados terminados de baja calidad (todos los pisos)
LotShape forma general de la propiedad GrLivArea pies cuadrados de superficie habitable sobre el nivel del suelo (suelo)
LandContour planitud de la propiedad BsmtFullBath baños completos del sótano
Utilities tipos de utilidades disponibles BsmtHalfBath medios baños del sótano
LotConfig configuración de las propiedades FullBath baños completos sobre el nivel del suelo
LandSlope pendiente de la propiedad HalfBath medios baños sobre el nivel del suelo
Neighborhood ubicaciones físicas dentro de los límites de la ciudad de Ames Bedroom número de dormitorios por encima del nivel del sótano
Condition1 proximidad a carretera principal o vía férrea GarageFinish acabado interior del garaje
Condition2 proximidad a la carretera principal o vía férrea (si hay una segunda presente) GarageCars tamaño del garaje en capacidad de automóviles
BldgType tipo de vivienda GarageArea tamaño del garaje en pies cuadrados
HouseStyle estilo de vivienda GarageQual calidad de garaje
OverallQual calidad general del material y del acabado GarageCond estado del garaje
OverallCond calificación de las condiciones generales PavedDrive entrada pavimentada
YearBuilt fecha de construcción original WoodDeckSF área de la cubierta de madera en pies cuadrados
YearRemodAdd fecha de remodelación OpenPorchSF área de porche abierto en pies cuadrados
RoofStyle tipo de techo EnclosedPorch área de porche cerrado en pies cuadrados
RoofMatl material del techo 3SsnPorch área de porche de tres estaciones en pies cuadrados
Exterior1st revestimiento exterior de la casa ScreenPorch área de porche de pantalla en pies cuadrados
Exterior2nd revestimiento exterior de la casa (si hay más de un material) PoolArea área de la piscina en pies cuadrados
MasVnrType tipo de chapa de mampostería PoolQC calidad de la piscina
MasVnrArea área de revestimiento de mampostería en pies cuadrados Fence calidad de la cerca
ExterQual calidad del material exterior MiscFeature miscelánea no incluida en otras categorías
ExterCond estado actual del material en el exterior MiscVal miscelánea referente al valor
Foundation tipo de cimentación MoSold mes de venta
BsmtQual altura del sótano YrSold año de venta
BsmtCond estado general del sótano SaleType tipo de venta
BsmtExposure paredes de los sótanos a nivel del jardín o de la salida SaleCondition condiciones de venta
BsmtFinType1 calidad del área terminada del sótano Kitchen número de cocinas
BsmtFinSF1 tipo 1: pies cuadrados terminados KitchenQual calidad de la cocina
BsmtFinType2 calidad de la segunda área terminada (si está presente) TotRmsAbvGrd total de habitaciones sobre rasante (no incluye baños)
BsmtFinSF2 tipo 2: pies cuadrados terminados Functional calificación de la funcionalidad del hogar
BsmtUnfSF pies cuadrados sin terminar del área del sótano Fireplaces número de chimeneas
TotalBsmtSF pies cuadrados totales del área del sótano FireplaceQu calidad de chimenea
GarageType ubicación del garaje GarageYrBlt año en que se construyó el garaje

Distribución de la variable objetivo

El objetivo será predecir el precio de una vivienda en función de sus características, por lo que SalesPrice será nuestra variable objetivo. En primer lugar comprobaremos cómo se distribuyen los valores de la objetivo.

# Objetivo: predecir el precio de una vivienda en función de sus características
house_complete |> 
  ggplot(aes(x = SalePrice)) +
  geom_density(alpha = .8, fill="#EB9891") +
  labs(title = "Distribución del precio de las viviendas", x = "Precio", y = NULL) +
  theme_minimal() +
  theme(text = element_text(face = "bold", size = 15), plot.title = element_text(hjust = 0.5), 
        axis.title.x = element_text(vjust=-0.5)) + 
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_y_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_vline(aes(xintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) +
  geom_vline(aes(xintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) + 
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted"))

En este caso, como la distribución parece bastante asimétrica, la mediana será la medida más representativa. Esta se encuentra en torno a los 160 000 dólares por vivienda.

Fase 2: Exploración de los datos

Problemas de codificación

Tras esta pequeña aproximación al dataset, comienza la primera fase de la metodología SEMMA para el depurado de nuestros datos. En este primer apartado observaremos a grosso modo si existen problemas de codificación en el dataset. Lo que más llama la atención con respecto a la codificación de las modalidades de ciertas variables es el gran número de ausentes que presenta el dataset.

ausentes <- 
  apply(house_complete, 2, function(x) sum(is.na(x)))

ausentes_tb <- 
  tibble(Variable = names(house_complete), Ausentes = ausentes) |> 
  filter(Ausentes > 0)
ausentes_tb
# A tibble: 20 × 2
   Variable    Ausentes
   <chr>          <int>
 1 MSZoning           4
 2 LotFrontage      486
 3 Alley           2721
 4 Utilities          2
 5 Exterior1st        1
 6 BsmtQual          81
 7 BsmtCond          82
 8 TotalBsmtSF        1
 9 Electrical         1
10 KitchenQual        1
11 Functional         2
12 FireplaceQu     1420
13 GarageYrBlt      159
14 GarageCars         1
15 GarageArea         1
16 GarageQual       159
17 PoolQC          2909
18 Fence           2348
19 SaleType           1
20 SalePrice       1459
ausentes_tb |> 
  filter(Ausentes > 10) |> 
  ggplot(aes(x = Variable, y = Ausentes)) +
  geom_col(position = "dodge", fill="#EB9891", color = "black") +
  labs(title = "> 10 valores ausentes por variable", x = "Variable", 
       y = "Cantidad de valores ausentes") +
  theme_minimal() +
  theme(text = element_text(face = "bold", size = 9), plot.title = element_text(hjust = 0.5),
        axis.title.x = element_text(vjust=-0.5), axis.title.y = element_text(vjust=2)) +
  scale_x_discrete(guide = guide_axis(n.dodge = 2))

Como se puede observar, 20 del total de variables presenta algún dato ausente. La mitad de esos 20 presenta más de 10 ausentes entre sus categorías. Ante esta cantidad de valores vacíos, debemos consultar primero si se tratan de categorías reales mal codificadas o ausentes reales. Según el documento explicativo, las variables Alley, BsmtQual, BsmtCond, FireplaceQu, GarageQual, PoolQC y Fence tienen codificadas como valores ausentes categorías que realmente no lo son. Vamos a modificar esta cuestión.

house <-
  house_complete |> 
  mutate(Alley = 
           ifelse(is.na(Alley), "None", Alley)) |>
  mutate(BsmtQual = 
           ifelse(is.na(BsmtQual), "None", BsmtQual)) |> 
  mutate(BsmtCond = 
           ifelse(is.na(BsmtCond), "None", BsmtCond)) |>
  mutate(FireplaceQu = 
           ifelse(is.na(FireplaceQu), "None", FireplaceQu)) |> 
  mutate(GarageQual = 
           ifelse(is.na(GarageQual), "None", GarageQual)) |>
  mutate(PoolQC = 
           ifelse(is.na(PoolQC), "None", PoolQC)) |>
  mutate(Fence = 
           ifelse(is.na(Fence), "None", Fence))

ausentes <- 
  apply(house, 2, function(x) sum(is.na(x)))

ausentes_tb <- 
  tibble(Variable = names(house), Ausentes = ausentes) |> 
  filter(Ausentes > 0)
ausentes_tb
# A tibble: 13 × 2
   Variable    Ausentes
   <chr>          <int>
 1 MSZoning           4
 2 LotFrontage      486
 3 Utilities          2
 4 Exterior1st        1
 5 TotalBsmtSF        1
 6 Electrical         1
 7 KitchenQual        1
 8 Functional         2
 9 GarageYrBlt      159
10 GarageCars         1
11 GarageArea         1
12 SaleType           1
13 SalePrice       1459

Ahora sí se ha reducido el número de variables con ausentes de 20 a 13. Al resto, le imputaremos lo que le corresponda (media, moda o mediana) en la fase pertinente.

Variables tipo texto, variables numéricas y factores

Tras ello, comprobaremos que todas las variables estén codificadas en su tipología correcta: debemos decidir si las variables de tipo texto son realmente variables cualitativas (factores).

house_complete |>  
  select(where(is.character)) |>
  glimpse()
Rows: 2,919
Columns: 28
$ MSZoning      <chr> "RL", "RL", "RL", "RL", "RL", "RL", "RL", "RL"…
$ Street        <chr> "Pave", "Pave", "Pave", "Pave", "Pave", "Pave"…
$ Alley         <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LandContour   <chr> "Lvl", "Lvl", "Lvl", "Lvl", "Lvl", "Lvl", "Lvl…
$ Utilities     <chr> "AllPub", "AllPub", "AllPub", "AllPub", "AllPu…
$ LandSlope     <chr> "Gtl", "Gtl", "Gtl", "Gtl", "Gtl", "Gtl", "Gtl…
$ Neighborhood  <chr> "CollgCr", "Veenker", "CollgCr", "Crawfor", "N…
$ BldgType      <chr> "1Fam", "1Fam", "1Fam", "1Fam", "1Fam", "1Fam"…
$ HouseStyle    <chr> "2Story", "1Story", "2Story", "2Story", "2Stor…
$ RoofStyle     <chr> "Gable", "Gable", "Gable", "Gable", "Gable", "…
$ RoofMatl      <chr> "CompShg", "CompShg", "CompShg", "CompShg", "C…
$ Exterior1st   <chr> "VinylSd", "MetalSd", "VinylSd", "Wd Sdng", "V…
$ ExterCond     <chr> "TA", "TA", "TA", "TA", "TA", "TA", "TA", "TA"…
$ BsmtQual      <chr> "Gd", "Gd", "Gd", "TA", "Gd", "Gd", "Ex", "Gd"…
$ BsmtCond      <chr> "TA", "TA", "TA", "Gd", "TA", "TA", "TA", "TA"…
$ Heating       <chr> "GasA", "GasA", "GasA", "GasA", "GasA", "GasA"…
$ HeatingQC     <chr> "Ex", "Ex", "Ex", "Gd", "Ex", "Ex", "Ex", "Ex"…
$ CentralAir    <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "…
$ Electrical    <chr> "SBrkr", "SBrkr", "SBrkr", "SBrkr", "SBrkr", "…
$ KitchenQual   <chr> "Gd", "TA", "Gd", "Gd", "Gd", "TA", "Gd", "TA"…
$ Functional    <chr> "Typ", "Typ", "Typ", "Typ", "Typ", "Typ", "Typ…
$ FireplaceQu   <chr> NA, "TA", "TA", "Gd", "TA", NA, "Gd", "TA", "T…
$ GarageQual    <chr> "TA", "TA", "TA", "TA", "TA", "TA", "TA", "TA"…
$ PavedDrive    <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "…
$ PoolQC        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ Fence         <chr> NA, NA, NA, NA, NA, "MnPrv", NA, NA, NA, NA, N…
$ SaleType      <chr> "WD", "WD", "WD", "WD", "WD", "WD", "WD", "WD"…
$ SaleCondition <chr> "Normal", "Normal", "Normal", "Abnorml", "Norm…

En el dataset encontramos varias variables divididas en modalidades del tipo Ex (excelente) > Gd (bueno) > TA (normal) > Fa (regular) > Po (pésimo). En concreto, se tratan de las variables ExterCond, BsmtQual, BsmtCond, HeatingQC, KitchenQual, FireplaceQu, GarageQual y PoolQC. Para ampliar el análisis de correlaciones posterior para con nuestra variable objetivo, lo que haremos será transformar todas estas variables cualitativas en ordinales. Como la mayoría sigue la misma estructura, crearemos un vector que asigne a cada modalidad un número del 0 al 5, siendo 0 igual a nada, y siendo 5 igual a excelente:

quality <- 
  c("None" = 0, "Po" = 1, "Fa" = 2, "TA" = 3, "Gd" = 4, "Ex" = 5)

Aplicamos el vector con revalue() a cada una de las variables:

house$ExterCond <-
  as.integer(revalue(house$ExterCond, quality))
house$BsmtQual <-
  as.integer(revalue(house$BsmtQual, quality))
house$BsmtCond <-
  as.integer(revalue(house$BsmtCond, quality))
house$HeatingQC <-
  as.integer(revalue(house$HeatingQC, quality))
house$KitchenQual <-
  as.integer(revalue(house$KitchenQual, quality))
house$FireplaceQu <-
  as.integer(revalue(house$FireplaceQu, quality))
house$GarageQual <-
  as.integer(revalue(house$GarageQual, quality))
house$PoolQC <-
  as.integer(revalue(house$PoolQC, quality))

Para el resto de variables no ordinales las convertimos a factor.

house <-
  house |> 
  mutate_if(~!is.numeric(.), as.factor)

Daremos un repaso también a las numéricas.

house_complete |>  
  select(where(is.numeric)) |>
  glimpse()
Rows: 2,919
Columns: 24
$ Id           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, …
$ MSSubClass   <dbl> 60, 20, 60, 70, 60, 50, 20, 60, 50, 190, 20, 60…
$ LotFrontage  <dbl> 65, 80, 68, 60, 84, 85, 75, NA, 51, 50, 70, 85,…
$ LotArea      <dbl> 8450, 9600, 11250, 9550, 14260, 14115, 10084, 1…
$ OverallCond  <dbl> 5, 8, 5, 5, 5, 5, 5, 6, 5, 6, 5, 5, 6, 5, 5, 8,…
$ YearBuilt    <dbl> 2003, 1976, 2001, 1915, 2000, 1993, 2004, 1973,…
$ YearRemodAdd <dbl> 2003, 1976, 2002, 1970, 2000, 1995, 2005, 1973,…
$ TotalBsmtSF  <dbl> 856, 1262, 920, 756, 1145, 796, 1686, 1107, 952…
$ `1stFlrSF`   <dbl> 856, 1262, 920, 961, 1145, 796, 1694, 1107, 102…
$ `2ndFlrSF`   <dbl> 854, 0, 866, 756, 1053, 566, 0, 983, 752, 0, 0,…
$ GrLivArea    <dbl> 1710, 1262, 1786, 1717, 2198, 1362, 1694, 2090,…
$ FullBath     <dbl> 2, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 3, 1, 2, 1, 1,…
$ HalfBath     <dbl> 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,…
$ BedroomAbvGr <dbl> 3, 3, 3, 3, 4, 1, 3, 3, 2, 2, 3, 4, 2, 3, 2, 2,…
$ TotRmsAbvGrd <dbl> 8, 6, 6, 7, 9, 5, 7, 7, 8, 5, 5, 11, 4, 7, 5, 5…
$ Fireplaces   <dbl> 0, 1, 1, 1, 1, 0, 1, 2, 2, 2, 0, 2, 0, 1, 1, 0,…
$ GarageYrBlt  <dbl> 2003, 1976, 2001, 1998, 2000, 1993, 2004, 1973,…
$ GarageCars   <dbl> 2, 2, 2, 3, 3, 2, 2, 2, 2, 1, 1, 3, 1, 3, 1, 2,…
$ GarageArea   <dbl> 548, 460, 608, 642, 836, 480, 636, 484, 468, 20…
$ OpenPorchSF  <dbl> 61, 0, 42, 35, 84, 30, 57, 204, 0, 4, 0, 21, 0,…
$ PoolArea     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ MoSold       <dbl> 2, 5, 9, 2, 12, 10, 8, 11, 4, 1, 2, 7, 9, 8, 5,…
$ YrSold       <dbl> 2008, 2007, 2008, 2006, 2008, 2009, 2007, 2009,…
$ SalePrice    <dbl> 208500, 181500, 223500, 140000, 250000, 143000,…

Nos llama la atención la variable MSSubClass, que identifica el tipo de vivienda objeto de venta. Aunque este codificada como numérica, realmente se trata de una variable cualitativa, por lo que vamos a convertirla.

class <- 
  c("20" = "1 story 1946+", "30" = "1 story 1945-", 
    "40" = "1 story fin attic", "45" = "1,5 story unf", 
    "50" = "1,5 story fin", "60" = "2 story 1946+", 
    "70" = "2 story 1945-", "75" = "2,5 story all ages", 
    "80" = "split/multi level", "85" = "split foyer", 
    "90" = "duplex all style/age", "120" = "1 story PUD 1946+", 
    "150" = "1,5 story PUD all", "160" = "2 story PUD 1946+", 
    "180" = "PUD multilevel", "190" = "2 family conversion")

house$MSSubClass <-
  as.factor(house$MSSubClass)
house$MSSubClass <-
  revalue(house$MSSubClass, class)

Ahora sí, todas las variables cualitativas y cuantitativas están bien codificadas.

Variables cuantitativas

Una vez asignado a cada variable su tipología correspondiente, pasaremos a analizar las variables cuantitativas del dataset. Se analizará ante todo cómo afecta cada variable a nuestra variable objetivo (SalePrice). Este análisis servirá, ante todo, para recategorizar las predictoras numéricas a la hora de crear la receta para el árbol en modo regresión.

Colinealidad

Antes de nada, al tener que analizar un dataset con tantas variables, comprobaremos los posibles problemas de colinealidad entre las predictoras numéricas con tal de eliminar las que repitan información. También, al tener una variable continua como objetivo, comprobaremos cuáles de las numéricas tienen una mayor correlación con ella con tal de mantenerlas y analizarlas en profundidad.

library(corrr)
cor_matrix <- 
  house |> select(where(is.numeric)) |> cor(use = "pairwise.complete.obs", method = "pearson")
library(corrplot)
cor_matrix |>
  corrplot(method = "number", tl.cex = 0.55, number.cex = 0.7, type = "lower")

Como se puede observar, existen bastantes variables con correlaciones superiores a 0.49. Vamos a echarles un vistazo:

# Seleccionar solo las entradas de la matriz de correlación con valores mayores a 0.49
filtered_matrix <- abs(cor_matrix)
filtered_matrix[cor_matrix <= 0.49] <- NA

# Creamos el gráfico de correlación
corrplot(filtered_matrix, method = "number", type = "lower",tl.cex = 0.55, number.cex = 0.7, na.label=" ")

En primer lugar, a partir del gráfico podemos observar como las variables YearsBuilt, YearRemodAdd, BsmtQual, TotalBsmtSF, 1stFlrSF, GrLivArea, FullBath, KitchenQual, TotRmsAbvGrd, FireplaceQu GarageCars y GarageArea son las que mantienen correlaciones más fuertes (superiores a 0.5 en valor absoluto) con nuestra variable objetivo SalePrice. Podríamos añadir también las variables Fireplaces y GarageYrBlt ya que mantienen correlaciones de 0.47 y 0.49, respectivamente.

En segundo lugar, analizaremos las correlaciones entre las propias predictoras con tal de deshacernos de información redundante. Nos desharemos de aquellas variables que mantengan correlaciones superiores a 0.7 con otras variables. De un par de variables muy correlacionadas, eliminaremos la que menor correlación mantenga con nuestra objetivo (SalePrice).

house <- 
  house |> 
  select(-c(GarageYrBlt, TotalBsmtSF, TotRmsAbvGrd, Fireplaces, GarageArea, PoolArea))

house |> 
  select(where(is.numeric)) |> 
  glimpse()
Rows: 2,919
Columns: 25
$ Id           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, …
$ LotFrontage  <dbl> 65, 80, 68, 60, 84, 85, 75, NA, 51, 50, 70, 85,…
$ LotArea      <dbl> 8450, 9600, 11250, 9550, 14260, 14115, 10084, 1…
$ OverallCond  <dbl> 5, 8, 5, 5, 5, 5, 5, 6, 5, 6, 5, 5, 6, 5, 5, 8,…
$ YearBuilt    <dbl> 2003, 1976, 2001, 1915, 2000, 1993, 2004, 1973,…
$ YearRemodAdd <dbl> 2003, 1976, 2002, 1970, 2000, 1995, 2005, 1973,…
$ ExterCond    <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,…
$ BsmtQual     <int> 4, 4, 4, 3, 4, 4, 5, 4, 3, 3, 3, 5, 3, 4, 3, 3,…
$ BsmtCond     <int> 3, 3, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,…
$ HeatingQC    <int> 5, 5, 5, 4, 5, 5, 5, 5, 4, 5, 5, 5, 3, 5, 3, 5,…
$ `1stFlrSF`   <dbl> 856, 1262, 920, 961, 1145, 796, 1694, 1107, 102…
$ `2ndFlrSF`   <dbl> 854, 0, 866, 756, 1053, 566, 0, 983, 752, 0, 0,…
$ GrLivArea    <dbl> 1710, 1262, 1786, 1717, 2198, 1362, 1694, 2090,…
$ FullBath     <dbl> 2, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 3, 1, 2, 1, 1,…
$ HalfBath     <dbl> 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,…
$ BedroomAbvGr <dbl> 3, 3, 3, 3, 4, 1, 3, 3, 2, 2, 3, 4, 2, 3, 2, 2,…
$ KitchenQual  <int> 4, 3, 4, 4, 4, 3, 4, 3, 3, 3, 3, 5, 3, 4, 3, 3,…
$ FireplaceQu  <int> 0, 3, 3, 4, 3, 0, 4, 3, 3, 3, 0, 4, 0, 4, 2, 0,…
$ GarageCars   <dbl> 2, 2, 2, 3, 3, 2, 2, 2, 2, 1, 1, 3, 1, 3, 1, 2,…
$ GarageQual   <int> 3, 3, 3, 3, 3, 3, 3, 3, 2, 4, 3, 3, 3, 3, 3, 3,…
$ OpenPorchSF  <dbl> 61, 0, 42, 35, 84, 30, 57, 204, 0, 4, 0, 21, 0,…
$ PoolQC       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ MoSold       <dbl> 2, 5, 9, 2, 12, 10, 8, 11, 4, 1, 2, 7, 9, 8, 5,…
$ YrSold       <dbl> 2008, 2007, 2008, 2006, 2008, 2009, 2007, 2009,…
$ SalePrice    <dbl> 208500, 181500, 223500, 140000, 250000, 143000,…

Variable Id

house |> 
  summarise(min_lead = min(Id), max_lead = max(Id))
  min_lead max_lead
1        1     2919

Como se puede observar, la variable Id es un id del número de viviendas que se incluyen en el dataset. Como no tiene interés alguno a fin de predecir nuestra variable objetivo, en la fase de modificación eliminaremos esta variable. La recuperaremos únicamente para cuando llegue el momento de subir las predicciones a Kaggle.

Variables relacionadas con medidas de áreas de la vivienda

sum1 <- house |> 
  summarise(Variable = "GrLivArea", min_lead = min(GrLivArea), max_lead = max(GrLivArea))

sum2 <- house |> 
  summarise(Variable = "`1stFlrSF`", min_lead = min(`1stFlrSF`), max_lead = max(`1stFlrSF`))

sum3 <- house |> 
  summarise(Variable = "`2ndFlrSF`", min_lead = min(`2ndFlrSF`), max_lead = max(`2ndFlrSF`))

sum4 <- house |> 
  summarise(Variable = "LotArea", min_lead = min(LotArea), max_lead = max(LotArea))

sum5 <- house |> 
  drop_na(LotFrontage) |> 
  summarise(Variable = "LotFrontage", min_lead = min(LotFrontage), max_lead = max(LotFrontage))

sum6 <- house |> 
  summarise(Variable = "OpenPorchSF", min_lead = min (OpenPorchSF), max_lead = max(OpenPorchSF))

rbind(sum1, sum2, sum3, sum4, sum5, sum6)
     Variable min_lead max_lead
1   GrLivArea      334     5642
2  `1stFlrSF`      334     5095
3  `2ndFlrSF`        0     2065
4     LotArea     1300   215245
5 LotFrontage       21      313
6 OpenPorchSF        0      742

Las seis variables miden áreas de ciertas zonas de la vivienda. Las variables 2ndFlrSF y OpenPorchSF toman el valor 0 cuando la vivienda no tiene segunda planta o porche. Graficaremos a continuación las seis variables para comprobar como se distribuyen en función de nuestra objetivo SalePrice.

a1 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = GrLivArea, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_text_repel(aes(label = ifelse(house$GrLivArea[!is.na(house$SalePrice)] > 4500, 
                                     rownames(house), ""))) +
  theme_minimal()

a2 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = `1stFlrSF`, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_text_repel(aes(label = ifelse(house$`1stFlrSF`[!is.na(house$SalePrice)] > 4500, 
                                     rownames(house), ""))) +
  theme_minimal()

a3 <- house[!is.na(house$SalePrice), ] |> 
  filter(`2ndFlrSF` > 0) |> 
  ggplot(aes(x = `2ndFlrSF`, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

a4 <- house[!is.na(house$SalePrice), ] |> 
  filter(LotArea < 100000) |> 
  ggplot(aes(x = LotArea, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

a5 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = LotFrontage, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_text_repel(aes(label = ifelse(house$LotFrontage[!is.na(house$SalePrice)] > 250, 
                                     rownames(house), ""))) +
  theme_minimal()

a6 <- house[!is.na(house$SalePrice), ] |> 
  filter(OpenPorchSF > 0) |> 
  ggplot(aes(x = OpenPorchSF, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

multiplot(a1, a2, a3, a4, a5, a6, cols = 2)

Como se puede observar, no todas las variables evolucionan de la misma manera al relacionarlas con la variable SalePrice.

Como ya tratamos los ausentes de estas variables (a excepción de LotFrontage, a los que probablemente imputemos la mediana o la eliminemos por su elevado número de ausentes), dejaremos estas seis variables tal y como están, a excepción de OpenPorchSF, la cual seguramente eliminemos.

Variables relacionadas con calidades y condiciones de la vivienda

sum1 <- house |> 
  summarise(Variable = "OverallCond", min_lead = min(OverallCond), max_lead = max(OverallCond))

sum2 <- house |> 
  summarise(Variable = "ExterCond", min_lead = min(ExterCond), max_lead = max(ExterCond))

sum3 <- house |> 
  summarise(Variable = "BsmtQual", min_lead = min(BsmtQual), max_lead = max(BsmtQual))

sum4 <- house |> 
  summarise(Variable = "BsmtCond", min_lead = min(BsmtCond), max_lead = max(BsmtCond))

sum5 <- house |> 
  drop_na(KitchenQual) |> 
  summarise(Variable = "KitchenQual", min_lead = min(KitchenQual), max_lead = max(KitchenQual))

sum6 <- house |> 
  summarise(Variable = "GarageQual", min_lead = min(GarageQual), max_lead = max(GarageQual))

sum7 <- house |> 
  summarise(Variable = "FireplaceQu", min_lead = min(FireplaceQu), max_lead = max(FireplaceQu))

sum8 <- house |> 
  summarise(Variable = "PoolQC", min_lead = min(PoolQC), max_lead = max(PoolQC))

sum9 <- house |> 
  summarise(Variable = "HeatingQC", min_lead = min(HeatingQC), max_lead = max(HeatingQC))

rbind(sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9)
     Variable min_lead max_lead
1 OverallCond        1        9
2   ExterCond        1        5
3    BsmtQual        0        5
4    BsmtCond        0        4
5 KitchenQual        2        5
6  GarageQual        0        5
7 FireplaceQu        0        5
8      PoolQC        0        5
9   HeatingQC        1        5

Las nueve variables miden las condiciones y calidades de ciertas zonas de la vivienda. Las variables BsmtQual, GarageQual, FireplaceQu y PoolQC toman el valor 0 cuando la vivienda no dispone de esas áreas. Graficaremos a continuación las nueve variables para comprobar como se distribuyen en función de nuestra objetivo SalePrice.

c1 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(OverallCond), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Overall Condition") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c2 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(ExterCond), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Exterior Condition") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c3 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(BsmtQual), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Basement Quality") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c4 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(BsmtCond), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Basement Condition") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c5 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(KitchenQual), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Kitchen Quality") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c6 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(GarageQual), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Garage Quality") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c7 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(FireplaceQu), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Fireplace Quality") + 
  scale_y_continuous(breaks= seq(0, 800000, by=100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c8 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(PoolQC), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Pool Quality") + 
  scale_y_continuous(breaks= seq(0, 800000, by=100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c9 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(HeatingQC), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Heating Quality") + 
  scale_y_continuous(breaks= seq(0, 800000, by=100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

multiplot(c1, c2, c3, c4, c5, c6, c7, c8, c9, cols = 3)

A priori, podríamos pensar que a mayor calidad en las zonas de la vivienda, mayor será el precio de la misma, pero no con todas parece cumplirse esta afirmación.

b1 <- ggplot(house, aes(x = factor(OverallCond))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="OverallCond") +
  theme_minimal()

b2 <- ggplot(house, aes(x = factor(ExterCond))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="ExterCond") +
  theme_minimal()

b3 <- ggplot(house, aes(x = factor(BsmtQual))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="BsmtQual") +
  theme_minimal()

b4 <- ggplot(house, aes(x = factor(BsmtCond))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="BsmtCond") +
  theme_minimal()

b5 <- ggplot(house, aes(x = factor(KitchenQual))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="KitchenQual") +
  theme_minimal()

b6 <- ggplot(house, aes(x = factor(GarageQual))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="GarageQual") +
  theme_minimal()

b7 <- ggplot(house, aes(x = factor(FireplaceQu))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="FireplaceQu") +
  theme_minimal()

b8 <- ggplot(house, aes(x = factor(PoolQC))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="PoolQC") +
  theme_minimal()

b9 <- ggplot(house, aes(x = factor(HeatingQC))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="HeatingQC") +
  theme_minimal()

multiplot(b1, b2, b3, b4, b5, b6, b7, b8, b9, cols = 3)

Por otro lado, si se atiende a la representación de cada categoría sobre el total de registros, en las variables GarageQual, BsmtCond o ExterCond la mayor parte las gobiernan la modalidad 3. Posteriormente, en la fase de la preparación de la receta, se podrían agrupar algunos de estos niveles con tal de romper la concentración de los datos o, por el contrario, eliminar directamente este tipo de variables.

De momento, mantendremos todas las variables a excepción de PoolQC, que será sustituida por una binaria.

Variables relacionadas con la edad de la vivienda

sum1 <- house |> 
  summarise(Variable = "YearRemodAdd", min_lead = min(YearRemodAdd), max_lead = max(YearRemodAdd))

sum2 <- house |> 
  summarise(Variable = "YearBuilt", min_lead = min(YearBuilt), max_lead = max(YearBuilt))

sum3 <- house |> 
  summarise(Variable = "MoSold", min_lead = min(MoSold), max_lead = max(MoSold))

sum4 <- house |> 
  summarise(Variable = "YrSold", min_lead = min(YrSold), max_lead = max(YrSold))

rbind(sum1, sum2, sum3, sum4)
      Variable min_lead max_lead
1 YearRemodAdd     1950     2010
2    YearBuilt     1872     2010
3       MoSold        1       12
4       YrSold     2006     2010

Las cuatro variables miden instantes temporales (años o meses) sobre distintas cuestiones relacionadas con la vida de la vivienda. La variable MoSold toma valores del 1 al 12 porque informa únicamente del mes del año en el que se vendió la vivienda. Graficaremos a continuación las cuatro variables para comprobar como se distribuyen en función de nuestra objetivo SalePrice.

t1 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = YearRemodAdd, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

t2 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = YearBuilt, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

t3 <- aggregate(SalePrice ~ YrSold, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$YrSold))) |> 
  ggplot(aes(x = YrSold, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#EB9891") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 25000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  coord_cartesian(ylim = c(0, 200000)) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Year Sold",y = "Sale Price") +
  theme_minimal()

t4 <- aggregate(SalePrice ~ MoSold, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$MoSold))) |> 
  ggplot(aes(x = MoSold, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#EB9891") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 800000, by = 25000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  coord_cartesian(ylim = c(0, 200000)) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Month Sold",y = "Sale Price") +
  theme_minimal()

multiplot(t1, t2, t3, t4, cols = 2)

Algunos comentarios a la vista de los gráficos:

En resumen, con estas cuatro variables crearemos otras dos nuevas para, posteriormente, eliminar la variable YearRemodAdd original. Además, las variables YrSold y MoSold podrían funcionar también como cualitativas y factores.

house$YrSold <-
  as.factor(house$YrSold)
house$MoSold <-
  as.factor(house$MoSold)

Variables relacionadas con los baños de la vivienda

sum1 <- house |> 
  summarise(Variable = "FullBath", min_lead = min(FullBath), max_lead = max(FullBath))

sum2 <- house |> 
  summarise(Variable = "HalfBath", min_lead = min(HalfBath), max_lead = max(HalfBath))

rbind(sum1, sum2)
  Variable min_lead max_lead
1 FullBath        0        4
2 HalfBath        0        2

Estas dos variables cuentan los baños totales y los medios baños de la vivienda. Estas toman el valor 0 cuando la vivienda no dispone de estos tipos de baño. Graficaremos a continuación las dos variables para comprobar como se distribuyen en función de nuestra objetivo SalePrice.

c1 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(FullBath), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Full Bath") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c2 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(HalfBath), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Half Bath") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

multiplot(c1, c2)

Como ninguna de estas dos variables parece muy potente, ni tampoco tiene una influencia determinante sobre nuestra objetivo, lo que se ha decidido es unirlas para conformar una nueva variable. De esta manera, la suma de FullBath y \(HalfBath * 0.5\) nos ofrecería información acerca del total de baños en la vivienda. La variable HalfBath (medios baños) se ponderaría por su significado real dividiendo sus valores y, por tanto, su importancia entre la mitad.

Variables restantes

sum1 <- house |> 
  summarise(Variable = "BedroomAbvGr", min_lead = min(BedroomAbvGr), max_lead = max(BedroomAbvGr))

sum2 <- house |> 
  drop_na(GarageCars) |> 
  summarise(Variable = "GarageCars", min_lead = min(GarageCars), max_lead = max(GarageCars))

rbind(sum1, sum2)
      Variable min_lead max_lead
1 BedroomAbvGr        0        8
2   GarageCars        0        5

Estas dos variables cuentan las habitaciones y los estacionamientos del garaje de la vivienda. Estas toman el valor 0 cuando la vivienda no dispone de estas áreas. Graficaremos a continuación las dos variables para comprobar como se distribuyen en función de nuestra objetivo SalePrice.

c1 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(BedroomAbvGr), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Bedrooms") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

c2 <- house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = factor(GarageCars), y = SalePrice)) +
  geom_boxplot(col = "#EB9891") + labs(x = "Parking lots") + 
  scale_y_continuous(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  theme_minimal()

multiplot(c1, c2)

Visualicemos también el peso de cada modalidad sobre el total de registros.

b1 <- ggplot(house, aes(x = factor(BedroomAbvGr))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="Bedroom") +
  theme_minimal()

b2 <- ggplot(house, aes(x = factor(GarageCars))) +
  geom_histogram(stat = "Count", fill = "#EB9891") + 
  labs(x="Parking lots") +
  theme_minimal()

multiplot(b1, b2)

A la luz de los gráficos, se puede observar como ambas variables siguen una correlación positiva para con nuestra variable objetivo (conclusión coherente, como con las anteriores variables). A partir de los valores 3-4, el número de registros decae para esas modalidades. Seguramente se opte en la fase de reagrupación por eliminar o incluir esas modalidades en la categoría contigua.

Variables cualitativas

A continuación analizaremos y agruparemos las variables cualitativas de cara a la creación de la receta del algoritmo KNN en modo regresión.

Variables relacionadas con las partes exteriores de la vivienda

house |> 
  select(where(is.factor)) |> 
  glimpse()
Rows: 2,919
Columns: 23
$ MSSubClass    <fct> "2 story 1946+", "1 story 1946+", "2 story 194…
$ MSZoning      <fct> RL, RL, RL, RL, RL, RL, RL, RL, RM, RL, RL, RL…
$ Street        <fct> Pave, Pave, Pave, Pave, Pave, Pave, Pave, Pave…
$ Alley         <fct> None, None, None, None, None, None, None, None…
$ LandContour   <fct> Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, L…
$ Utilities     <fct> AllPub, AllPub, AllPub, AllPub, AllPub, AllPub…
$ LandSlope     <fct> Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, G…
$ Neighborhood  <fct> CollgCr, Veenker, CollgCr, Crawfor, NoRidge, M…
$ BldgType      <fct> 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, 1Fam, 1Fam…
$ HouseStyle    <fct> 2Story, 1Story, 2Story, 2Story, 2Story, 1.5Fin…
$ RoofStyle     <fct> Gable, Gable, Gable, Gable, Gable, Gable, Gabl…
$ RoofMatl      <fct> CompShg, CompShg, CompShg, CompShg, CompShg, C…
$ Exterior1st   <fct> VinylSd, MetalSd, VinylSd, Wd Sdng, VinylSd, V…
$ Heating       <fct> GasA, GasA, GasA, GasA, GasA, GasA, GasA, GasA…
$ CentralAir    <fct> Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y…
$ Electrical    <fct> SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SBrk…
$ Functional    <fct> Typ, Typ, Typ, Typ, Typ, Typ, Typ, Typ, Min1, …
$ PavedDrive    <fct> Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y…
$ Fence         <fct> None, None, None, None, None, MnPrv, None, Non…
$ MoSold        <fct> 2, 5, 9, 2, 12, 10, 8, 11, 4, 1, 2, 7, 9, 8, 5…
$ YrSold        <fct> 2008, 2007, 2008, 2006, 2008, 2009, 2007, 2009…
$ SaleType      <fct> WD, WD, WD, WD, WD, WD, WD, WD, WD, WD, WD, Ne…
$ SaleCondition <fct> Normal, Normal, Normal, Abnorml, Normal, Norma…
b1 <- aggregate(SalePrice ~ Street, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Street))) |> 
  ggplot(aes(x = Street, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Street",y = "Sale Price") +
  theme_minimal()

b2 <- aggregate(SalePrice ~ Alley, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Alley))) |> 
  ggplot(aes(x = Alley, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Alley",y = "Sale Price") +
  theme_minimal()

b3 <- aggregate(SalePrice ~ LandContour, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$LandContour))) |> 
  ggplot(aes(x = LandContour, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "LandContour",y = "Sale Price") +
  theme_minimal()

b4 <- aggregate(SalePrice ~ LandSlope, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$LandSlope))) |> 
  ggplot(aes(x = LandSlope, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "LandSlope",y = "Sale Price") +
  theme_minimal()

b5 <- aggregate(SalePrice ~ RoofStyle, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$RoofStyle))) |> 
  ggplot(aes(x = RoofStyle, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "RoofStyle",y = "Sale Price") +
  scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  theme_minimal()

b6 <- aggregate(SalePrice ~ RoofMatl, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$RoofMatl))) |> 
  ggplot(aes(x = RoofMatl, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.7)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "RoofMatl",y = "Sale Price") +
  coord_flip() +
  theme_minimal()

b7 <- aggregate(SalePrice ~ Exterior1st, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Exterior1st))) |> 
  ggplot(aes(x = Exterior1st, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.5)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Exterior1st",y = "Sale Price")   +
  coord_flip() +
  theme_minimal()

b8 <- aggregate(SalePrice ~ PavedDrive, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$PavedDrive))) |> 
  ggplot(aes(x = PavedDrive, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "PavedDrive",y = "Sale Price") +
  theme_minimal()

b9 <- aggregate(SalePrice ~ Fence, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Fence))) |> 
  ggplot(aes(x = Fence, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Fence",y = "Sale Price") +
  theme_minimal()

multiplot(b1, b2, b3, b4, b5, b6, b7, b8, b9, cols = 3)

A la luz de los gráficos, todas las variables presentan cierta influencia sobre nuestra variable objetivo, aunque muchas de sus categorías son muy minoritarias. Resulta necesario reagrupar la mayoría de variables. Se agruparan en función de su influencia sobre nuestra objetivo.

Variables relacionadas con la zona en la que se sitúa la vivienda

b1 <- aggregate(SalePrice ~ MSSubClass, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$MSSubClass))) |> 
  ggplot(aes(x = MSSubClass, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.7)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "MSSubClass",y = "Sale Price")  +
  coord_flip() +
  theme_minimal()

b2 <- aggregate(SalePrice ~ MSZoning, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$MSZoning))) |> 
  ggplot(aes(x = MSZoning, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "MSZoning",y = "Sale Price") +
  theme_minimal()

b3 <- aggregate(SalePrice ~ Neighborhood, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Neighborhood))) |> 
  ggplot(aes(x = Neighborhood, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.7)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Neighborhood",y = "Sale Price") +
  coord_flip() +
  theme_minimal()

b4 <- aggregate(SalePrice ~ BldgType, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$BldgType))) |> 
  ggplot(aes(x = BldgType, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "BldgType",y = "Sale Price") +
  theme_minimal()

b5 <- aggregate(SalePrice ~ HouseStyle, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$HouseStyle))) |> 
  ggplot(aes(x = HouseStyle, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "HouseStyle",y = "Sale Price") +
  theme_minimal()

multiplot(b1, b2, b3, b4, b5, cols = 2)

Variables relacionadas con los servicios y utilidades de la vivienda

b1 <- aggregate(SalePrice ~ Utilities, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Utilities))) |> 
  ggplot(aes(x = Utilities, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Utilities",y = "Sale Price") +
  theme_minimal()

b2 <- aggregate(SalePrice ~ Heating, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Heating))) |> 
  ggplot(aes(x = Heating, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Heating",y = "Sale Price") +
  theme_minimal()

b3 <- aggregate(SalePrice ~ CentralAir, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$CentralAir))) |> 
  ggplot(aes(x = CentralAir, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "CentralAir",y = "Sale Price") +
  theme_minimal()

b4 <- aggregate(SalePrice ~ Electrical, house[!is.na(house$SalePrice), ], median) |>
  ggplot(aes(x = Electrical, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Electrical",y = "Sale Price") +
  theme_minimal()

b5 <- aggregate(SalePrice ~ Functional, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$Functional))) |> 
  ggplot(aes(x = Functional, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "Functional",y = "Sale Price") +
  theme_minimal()

multiplot(b1, b2, b3, b4, b5, cols = 2)

Variables relacionadas con la venta de la vivienda

b1 <- aggregate(SalePrice ~ SaleType, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$SaleType))) |> 
  ggplot(aes(x = SaleType, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "SaleType",y = "Sale Price") +
  theme_minimal()

b2 <- aggregate(SalePrice ~ SaleCondition, house[!is.na(house$SalePrice), ], median) |>
  mutate(n = pull(count(house[!is.na(house$SalePrice), ]$SaleCondition))) |> 
  ggplot(aes(x = SaleCondition, y = SalePrice)) +
  geom_bar(stat = "identity", fill= "#56BCC2") +
  geom_label(aes(label = n, y = 0.1)) +
  scale_y_continuous(breaks= seq(0, 700000, by = 50000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_hline(aes(yintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) + 
  geom_hline(aes(yintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) +
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted")) +
  labs(x = "SaleCondition",y = "Sale Price") +
  theme_minimal()

multiplot(b1, b2)

Fase 1 y 3: Muestreo y modificación de los datos

Tras la fase de exploración de los datos, continuaremos con las fases de muestreo y modificación de los datos. Dado que nuestro dataset contiene tan solo 2919 observaciones (1460 para train y 1459 para test), no será necesario realizar muestreo (nos quedaríamos prácticamente sin filas si lo hacemos). Por otro lado, iremos fijando semilla a fin de poder sacar conclusiones de los resultados de nuestros distintos modelos.

En segundo lugar, para la fase de modificación de los datos, consideraremos dos apartados principales. Uno primero en donde se ejecutarán las modificaciones estructurales que afecten a toda las base de datos (transformar variables a factores, problemas de codificación o de rango, variables que no aportan, creación de variables en general, etc.), y uno segundo en donde se llevarán a cabo aquellas modificaciones que afecten a cada algoritmo en concreto a modo de receta (normalización para la métrica, recategorización, tratamiento de outliers/ausentes, dummyficación, etc.).

Modificaciones estructurales sobre conjuntos test y train/validación

En este caso, como las particiones ya están realizadas, vamos a aplicar al conjunto completo (train+ test) las modificaciones estructurales para que tengan la misma forma.

# Aplicamos las modificaciones estructurales al dataset completo (train + test)
## Tratamiento de ausentes y transformación de variables texto a factor
house_complete <- 
  house_complete |>   
  mutate(Alley = 
           ifelse(is.na(Alley), "None", Alley)) |>
  mutate(BsmtQual = 
           ifelse(is.na(BsmtQual), "None", BsmtQual)) |> 
  mutate(BsmtCond = 
           ifelse(is.na(BsmtCond), "None", BsmtCond)) |>
  mutate(FireplaceQu = 
           ifelse(is.na(FireplaceQu), "None", FireplaceQu)) |> 
  mutate(GarageQual = 
           ifelse(is.na(GarageQual), "None", GarageQual)) |>
  mutate(PoolQC = 
           ifelse(is.na(PoolQC), "None", PoolQC)) |>
  mutate(Fence = 
           ifelse(is.na(Fence), "None", Fence))

## Modificaciones en la tipología de las variables
house_complete$ExterCond <-
  as.integer(revalue(house_complete$ExterCond, quality))
house_complete$BsmtQual <-
  as.integer(revalue(house_complete$BsmtQual, quality))
house_complete$BsmtCond <-
  as.integer(revalue(house_complete$BsmtCond, quality))
house_complete$HeatingQC <-
  as.integer(revalue(house_complete$HeatingQC, quality))
house_complete$KitchenQual <-
  as.integer(revalue(house_complete$KitchenQual, quality))
house_complete$FireplaceQu <-
  as.integer(revalue(house_complete$FireplaceQu, quality))
house_complete$GarageQual <-
  as.integer(revalue(house_complete$GarageQual, quality))
house_complete$PoolQC <-
  as.integer(revalue(house_complete$PoolQC, quality))
house_complete$KitchenQual[is.na(house_complete$KitchenQual)] <- 3
house_complete$GarageCars[is.na(house_complete$GarageCars)] <- 2

## Creación de nuevas variables
house_complete$RemodAdd <- # variable binaria ¿reformada?
  ifelse(house_complete$YearBuilt == house_complete$YearRemodAdd, 0, 1)

house_complete$New <- # variable binaria ¿es nueva?
  ifelse(house_complete$YrSold == house_complete$YearBuilt, 1, 0)

house_complete$TotBath <- # variable baño total
  house_complete$FullBath + (house_complete$HalfBath * 0.5)
house_complete <- 
  house_complete |> 
  select(-c(FullBath, HalfBath))

house_complete$Age <- # variable edad de la vivienda
  house_complete$YrSold-house_complete$YearRemodAdd
house_complete <- 
  house_complete |> 
  select(-c(YearRemodAdd))

## Transformación del resto a factor
house_complete$MSSubClass <-
  as.factor(house_complete$MSSubClass)
house_complete$YrSold <-
  as.factor(house_complete$YrSold)
house_complete$MoSold <-
  as.factor(house_complete$MoSold)

house_complete <-
  house_complete |>
  mutate_if(~!is.numeric(.), as.factor)

## Eliminación de variables por problemas de colinealidad
house_complete <- 
  house_complete |> 
  select(-c(GarageYrBlt, TotalBsmtSF, TotRmsAbvGrd, Fireplaces, GarageArea, PoolArea))

## Eliminación de variables por insustanciales
house_complete <- 
  house_complete |> 
  select(-c(Id, OpenPorchSF, Street, RoofMatl, Utilities, Heating, BedroomAbvGr, PoolQC))

División de particiones

Conjuntos train y test

Tras homogeneizar los subconjuntos de train y test, los volvemos a subdividir tal y como se encontraban en origen.

house_train <- house_complete[1:1460,]
house_test <- house_complete[1461:2919,]

Modificaciones dentro de la receta

Receta para el algoritmo KNN en modo regresión

Aplicación de roles

Tras las particiones, definimos la receta, indicándole el conjunto donde tenemos validación y train, y enfrentamos nuestra variable objetivo SalePrice con todas las demás. Después, asignamos posibles roles, sujetos a modificación, que nos permitan diferenciar acciones entre las variables.

# Receta
house_rec_knn <-
  # Fórmula y datos
  recipe(data = house_train, SalePrice ~ .)|>
  # Roles
  add_role(where(is.factor), 
           new_role = "cualitativa") |> 
  add_role(where(is.numeric), 
           new_role = "cuantitativa") |> 
  add_role(RemodAdd, New, 
           new_role = "binaria") |> 
  add_role(LotFrontage, LotArea, `1stFlrSF`, `2ndFlrSF`, GrLivArea, Age, 
           new_role = "mediana") |> 
  add_role(all_nominal_predictors(), 
           new_role = "moda") |> 
  add_role(ExterCond, BsmtCond, HeatingQC, GarageQual, KitchenQual,
           new_role = "imputar mediana") |> 
  add_role(OverallCond, BsmtQual, FireplaceQu, TotBath, YearBuilt, GarageCars,
           new_role = "imputar media")

Reagrupación de las variables cualitativas

En este apartado, reagrupamos las variables cualitativas con step_mutate con lo definido en la fase de exploración. Para el algoritmo de clasificación KNN, cuanto más agrupadas estén este tipo de categorías, mejor.

house_rec_knn <- 
  house_rec_knn |> 
  step_mutate(Alley = 
                fct_collapse(Alley, 
                             Grvl_Pave = c("Grvl", "Pave"))) |> 
  step_mutate(LandSlope = 
                fct_collapse(LandSlope, 
                             Mod_Sev = c("Mod", "Sev"))) |>  
  step_mutate(LandContour = 
                fct_collapse(LandContour, 
                             Bnk_HLS_Low = c("Bnk", "HLS", "Low"))) |> 
  step_mutate(RoofStyle = 
                fct_collapse(RoofStyle, 
                             Hip = c("Flat", "Shed", "Mansard", "Hip"),
                             Gam_Gable = c("Gambrel", "Gable"))) |> 
  step_mutate(Exterior1st = 
                fct_collapse(Exterior1st, 
                            `<150k` = c("AsphShn", "BrkComm", "BrkFace",
                                        "CemntBd", "HdBoard", "MetalSd",
                                        "Stucco", "Wd Sdng", "WdShing",
                                        "AsbShng"),
                             `>150k` = c("CBlock", "ImStucc", "Plywood",
                                         "Stone", "VinylSd"))) |> 
  step_mutate(PavedDrive = 
                fct_collapse(PavedDrive, 
                             N_P = c("N", "P"))) |>
  step_mutate(Fence = 
                fct_collapse(Fence, 
                             MnPrv = c("GdWo", "MnWw"),
                             None = c("GdPrv"))) |> 
  step_mutate(MSSubClass = 
                fct_collapse(MSSubClass, 
                             `<150k` = c("30", "40", "45", "50",
                                         "80", "85", "120", "160",
                                         "180", "190"),
                             `>150k` = c("20", "60", "70", "75",
                                         "90", "150"))) |>
  step_mutate(MSZoning = 
                fct_collapse(MSZoning, 
                             RM = c("C (all)", "RH"),
                             RL = c("FV"))) |> 
  step_mutate(Neighborhood = 
                fct_collapse(Neighborhood, 
                             `<125k` = c("BrDale", "BrkSide", "Edwards",
                                         "IDOTRR", "MeadowV", "OldTown"),
                             `125k-150k` = c("Blmngtn", "NAmes", "NPkVill",
                                             "SWISU", "Sawyer"),
                             `150k-200k` = c("Blueste", "CollgCr", "Gilbert",
                                             "Mitchel", "NWAmes", "SawyerW"),
                             `>200k` = c("ClearCr", "Crawfor", "NoRidge",
                                         "NridgHt", "Somerst", "StoneBr",
                                         "Timber", "Veenker"))) |> 
  step_mutate(BldgType = 
                fct_collapse(BldgType, 
                             Fam_TwnhsE = c("1Fam", "TwnhsE"),
                             fmCon_Dup_Twnhs = c("2fmCon", "Duplex", "Twnhs"))) |>  
  step_mutate(HouseStyle = 
                fct_collapse(HouseStyle, 
                             Story_Fin_SLvl = c("2Story", "2.5Fin", "SLvl"),
                             fmCon_Dup_Twnhs = c("1.5Fin", "1.5Unf", "1Story", 
                                                 "2.5Unf", "SFoyer"))) |> 
  step_mutate(Electrical = 
                fct_collapse(Electrical, 
                             Other = c("FuseA", "FuseF", "FuseP", "Mix"))) |>
  step_mutate(Functional = 
                fct_collapse(Functional, 
                             Other = c("Maj1", "Maj2", "Min1", "Min2", 
                                       "Mod", "Sev"))) |>
  step_mutate(SaleType = 
                fct_collapse(SaleType, 
                             Con_CWD_New = c("Con", "CWD", "New"),
                             WD = c("COD", "ConLD", "ConLI", "ConLw",
                                       "Oth"))) |> 
  step_mutate(SaleCondition = 
                fct_collapse(SaleCondition, 
                             Normal_Others = c("Abnorml", "AdjLand", "Alloca",
                                               "Family"))) |>
  step_mutate(MoSold = 
                fct_collapse(MoSold, 
                             Jan_Ap_May_Oct = c("1", "4", "5", "10"),
                             Mar_Jun_Jul = c("3", "6", "7"),
                             Feb_Aug_Sep_Nov_Dec = c("2", "8", "9",
                                                     "11", "12"))) |> 
  step_mutate(YrSold = 
                fct_collapse(YrSold, 
                             `2006-2009` = c("2006", "2007", "2008", "2009")))

Recategorización de las variables cuantitativas

En el algoritmo de clasificación KNN, la recategorización de las variables cuantitativas no es estrictamente necesaria. Por tanto, en este receta, del total de variables cuantitativas no vamos a modificar ninguna.

Outliers

box1 <- 
  ggplot(house_complete, aes(SalePrice, Age)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box2 <- 
  ggplot(house_complete, aes(SalePrice, LotFrontage)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box3 <- 
  ggplot(house_complete, aes(SalePrice, LotArea)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box4 <- 
  ggplot(house_complete, aes(SalePrice, OverallCond)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box5 <- 
  ggplot(house_complete, aes(SalePrice, YearBuilt)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box6 <- 
  ggplot(house_complete, aes(SalePrice, ExterCond)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box7 <- 
  ggplot(house_complete, aes(SalePrice, BsmtQual)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box8 <- 
  ggplot(house_complete, aes(SalePrice, BsmtCond)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box9 <- 
  ggplot(house_complete, aes(SalePrice, HeatingQC)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box10 <- 
  ggplot(house_complete, aes(SalePrice, `1stFlrSF`)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box11 <- 
  ggplot(house_complete, aes(SalePrice, `2ndFlrSF`)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box12 <- 
  ggplot(house_complete, aes(SalePrice, GrLivArea)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box13 <- 
  ggplot(house_complete, aes(SalePrice, KitchenQual)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box14 <- 
  ggplot(house_complete, aes(SalePrice, FireplaceQu)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box15 <- 
  ggplot(house_complete, aes(SalePrice, GarageCars)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box16 <- 
  ggplot(house_complete, aes(SalePrice, GarageQual)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()
box17 <- 
  ggplot(house_complete, aes(SalePrice, TotBath)) +
  geom_boxplot() +
  theme_minimal() +
  coord_flip()

multiplot(box1, box2, box3, box4, box5, box6, box7, box8, box9, box10, cols = 2) 
multiplot(box11, box12, box13, box14, box15, box16, box17, cols = 2)

Si observamos estos gráficos de cajas y bigotes, todas nuestras variables cuantitativas continuas son asimétricas, por lo que se detectarán los outliers y se imputarán los ausentes por la mediana. Para el caso de las semi-cuantitativas (KitchenQual, por ejemplo) se imputarán directamente los ausentes por uno de los dos métodos en función de la simetría o asimetría de la variable. Para el resto de variables cualitativas, imputamos directamente por la moda.

house_rec_knn <-
  house_rec_knn |> 
  # Detección de outliers por la mediana y por la media
  step_mutate(across(has_role("mediana"), function(x) { ifelse(abs(scores(x, type = "mad")) > 3 & !is.na(x), NA, x) })) |>
  # Imputación de ausentes por la mediana, la media y la moda
  step_impute_median(has_role("mediana")) |>
  step_impute_median(has_role("imputar_mediana")) |> 
  step_impute_mean(has_role("imputar_media")) |>
  step_impute_mode(has_role("moda"))

Filtro de correlación

Aplicamos el filtro de correlación a nuestras variables numéricas.

house_rec_knn <-
  house_rec_knn |> 
  step_corr(has_role("cuantitativa"), threshold = 0.9)

Normalizar por rango

En el caso del algoritmo KNN, necesitaremos normalizar nuestras variables por rango para que todas tengan el mismo peso, entre 0 y 1.

house_rec_knn <-
  house_rec_knn |> 
  step_range(all_numeric_predictors(), min = 0, max = 1) |> 
  step_other(all_nominal_predictors(), threshold = 0.02)

Variables dummy

Como esta receta es para el modelo KNN, debemos dummyficar nuestras variables cualitativas. Para ello, tomamos todas las nominales, menos nuestra variable objetivo.

house_rec_knn <-
  house_rec_knn |>
  step_dummy(all_nominal(), -all_outcomes())

Filtro de cero varianza

house_rec_knn <-
  house_rec_knn |>
  step_zv(all_predictors())

Horneado

Por último, horneamos nuestra receta para comprobar que todas nuestras nuevas variables recategorizadas se hayan creado correctamente.

bake(house_rec_knn |>  prep(), new_data = NULL)
# A tibble: 1,460 × 43
   LotFrontage LotArea OverallCond YearBuilt ExterCond BsmtQual
         <dbl>   <dbl>       <dbl>     <dbl>     <dbl>    <dbl>
 1      0.151   0.0334       0.5       0.949       0.5      0.8
 2      0.202   0.0388       0.875     0.754       0.5      0.8
 3      0.161   0.0465       0.5       0.935       0.5      0.8
 4      0.134   0.0386       0.5       0.312       0.5      0.6
 5      0.216   0.0606       0.5       0.928       0.5      0.8
 6      0.219   0.0599       0.5       0.877       0.5      0.8
 7      0.185   0.0411       0.5       0.957       0.5      1  
 8      0.164   0.0425       0.625     0.732       0.5      0.8
 9      0.103   0.0225       0.5       0.428       0.5      0.6
10      0.0993  0.0286       0.625     0.486       0.5      0.6
# … with 1,450 more rows, and 37 more variables: BsmtCond <dbl>,
#   HeatingQC <dbl>, `1stFlrSF` <dbl>, `2ndFlrSF` <dbl>,
#   GrLivArea <dbl>, KitchenQual <dbl>, FireplaceQu <dbl>,
#   GarageCars <dbl>, GarageQual <dbl>, RemodAdd <dbl>, New <dbl>,
#   TotBath <dbl>, Age <dbl>, SalePrice <dbl>,
#   MSSubClass_X.150k <dbl>, MSZoning_RL <dbl>, Alley_None <dbl>,
#   LandContour_Lvl <dbl>, LandSlope_Mod_Sev <dbl>, …

Receta para el método árboles de clasificación en modo regresión

Aplicación de roles

En este segundo apartado repetiremos la receta, pero esta vez adaptada a los árboles de clasificación.

# Receta
house_rec_arbol <-
  # Fórmula y datos
  recipe(data = house_train, SalePrice ~ .)|>
  # Roles
  add_role(where(is.factor), new_role = "cualitativa") |> 
  add_role(where(is.numeric), new_role = "cuantitativa") |> 
  add_role(RemodAdd, New, new_role = "binaria") |> 
  add_role(all_nominal_predictors(), new_role = "moda")

Reagrupación de las variables cualitativas

En este apartado, reagrupamos las variables cualitativas con step_mutate. Para el método de árboles de clasificación, no es necesario agrupar en exceso este tipo de variables, por lo que únicamente agruparemos las variables Neighborhood y Exterior1st (las que más modalidades presentaban).

house_rec_arbol <- 
  house_rec_arbol |> 
  step_mutate(Neighborhood = 
                fct_collapse(Neighborhood, 
                             `<125k` = c("BrDale", "BrkSide", "Edwards",
                                         "IDOTRR", "MeadowV", "OldTown"),
                             `125k-150k` = c("Blmngtn", "NAmes", "NPkVill",
                                             "SWISU", "Sawyer"),
                             `150k-200k` = c("Blueste", "CollgCr", "Gilbert",
                                             "Mitchel", "NWAmes", "SawyerW"),
                             `>200k` = c("ClearCr", "Crawfor", "NoRidge",
                                         "NridgHt", "Somerst", "StoneBr",
                                         "Timber", "Veenker"))) |> 
  step_mutate(Exterior1st = 
                fct_collapse(Exterior1st, 
                            `<150k` = c("AsphShn", "BrkComm", "BrkFace",
                                        "CemntBd", "HdBoard", "MetalSd",
                                        "Stucco", "Wd Sdng", "WdShing",
                                        "AsbShng"),
                             `>150k` = c("CBlock", "ImStucc", "Plywood",
                                         "Stone", "VinylSd")))

Recategorización de las variables cuantitativas

Como esta receta es para el método de árboles de clasificación, tenemos que recategorizar todas las variables cuantitativas.

house_rec_arbol <- 
  house_rec_arbol |> 
    step_mutate(GrLivArea =
                cut(GrLivArea,
                    breaks = c(0, 1000, 2000, 3000, Inf),
                    labels = c("<1000", "1000-2000", "2000-3000", ">3000"))) |> 
    step_mutate(`1stFlrSF` =
                cut(`1stFlrSF`,
                    breaks = c(0, 1000, 2000, Inf),
                    labels = c("<1000", "1000-2000", ">2000"))) |>
    step_mutate(`2ndFlrSF` =
                cut(`2ndFlrSF`,
                    breaks = c(0, 500, 1000, 1500, Inf),
                    labels = c("<500", "500-1000", "1000-1500", ">1500"))) |> 
    step_mutate(LotArea =
                cut(LotArea,
                    breaks = c(0, 10000, 20000, 30000, Inf),
                    labels = c("<10000", "10000-20000", "20000-30000", ">30000"))) |> 
    step_mutate(LotFrontage =
                cut(LotFrontage,
                    breaks = c(0, 50, 100, 150, Inf),
                    labels = c("<50", "50-100", "100-150", ">150"))) |> 
    step_mutate(OverallCond =
                cut(OverallCond,
                   breaks = c(-Inf, 3, 4, 5, 6, 7, Inf),
                   labels = c("1-3", "4", "5", "6", "7", "8-9"))) |> 
    step_mutate(BsmtCond =
                cut(BsmtCond,
                   breaks = c(-Inf, 1, 2, 3, Inf),
                   labels = c("0-1", "2", "3", "4"))) |> 
    step_mutate(FireplaceQu =
               cut(FireplaceQu,
                   breaks = c(-Inf, 2, 4, Inf),
                   labels = c("0-2", "3-4", "5"))) |> 
    step_mutate(ExterCond =
               cut(ExterCond,
                  breaks = c(-Inf, 2, 3, Inf),
                  labels = c("0-2", "3", "4-5"))) |> 
    step_mutate(KitchenQual =
               cut(KitchenQual,
                  breaks = c(-Inf, 3, 4, Inf),
                  labels = c("2-3", "4", "5"))) |> 
    step_mutate(BsmtQual =
               cut(BsmtQual,
                  breaks = c(-Inf, 2, 3, 4, Inf),
                  labels = c("0&2", "3", "4", "5"))) |> 
    step_mutate(GarageQual =
               cut(GarageQual,
                  breaks = c(-Inf, 1, 2, Inf),
                  labels = c("0-1", "2", "3-5"))) |> 
    step_mutate(HeatingQC =
               cut(HeatingQC,
                  breaks = c(-Inf, 3, 4, Inf),
                  labels = c("1-3", "4", "5"))) |> 
    step_mutate(YearBuilt =
               cut(YearBuilt,
                breaks = c(-Inf, 1920, 1960, Inf),
                labels = c("<1920", "1920-1960", ">1960"))) |>
    step_mutate(GarageCars =
               cut(GarageCars,
                  breaks = c(-Inf, 1, 2, 3, Inf),
                  labels = c("0-1", "2", "3", "4"))) |>
    step_mutate(Age =
               cut(Age,
                  breaks = c(-Inf, 10, 20, 30, 40, 50, Inf),
                  labels = c("0-10", "10-20", "20-30", "30-40", "40-50", "+50"))) |>
    step_mutate(TotBath =
               cut(TotBath,
                  breaks = c(-Inf, 1, 1.5, 2, 2.5, Inf),
                  labels = c("0-1", "1.5", "2", "2.5", "+3"))) |>
    step_mutate(RemodAdd =
               cut(RemodAdd,
                  breaks = c(-Inf, 0, Inf),
                  labels = c("No", "Yes"))) |>
    step_mutate(New =
               cut(New,
                  breaks = c(-Inf, 0, Inf),
                  labels = c("No", "Yes"))) |>
  step_mutate(across(where(is.character), as_factor))

Outliers

multiplot(box1, box2, box3, box4, box5, box6, box7, box8, box9, box10, cols = 2) 
multiplot(box11, box12, box13, box14, box15, box16, box17, cols = 2)

En esta receta se ha optado por recategorizar y reagrupar todas las variables, incluso las cuantitativas. De esta manera, no hay outliers que detectar, pues todas las nuevas categorías se han creado manualmente a partir de la agrupación de las antiguas modalidades. Aún así, imputamos a los posibles ausentes la moda al haber convertido todas las variables en cualitativas.

house_rec_arbol <-
  house_rec_arbol |> 
  step_impute_mode(all_predictors())

Filtro de correlación

Como todas nuestras variables están ahora recategorizadas y convertidas a factor, no existen variables estrictamente numéricas a las que aplicar el filtro de correlación.

Normalizar por rango

En árboles, no será necesario normalizar por rango.

Variables dummy

En árboles, tampoco será necesario dummyficar.

Filtro de cero varianza

house_rec_arbol <-
  house_rec_arbol |>
  step_zv(all_predictors())

Horneado

Por último, horneamos nuestra receta para comprobar que todas nuestras nuevas variables recategorizadas se hayan creado correctamente.

bake(house_rec_arbol |>  prep(), new_data = NULL)
# A tibble: 1,460 × 39
   MSSubClass MSZoning LotFrontage LotArea Alley LandContour LandSlope
   <fct>      <fct>    <fct>       <fct>   <fct> <fct>       <fct>    
 1 60         RL       50-100      <10000  None  Lvl         Gtl      
 2 20         RL       50-100      <10000  None  Lvl         Gtl      
 3 60         RL       50-100      10000-… None  Lvl         Gtl      
 4 70         RL       50-100      <10000  None  Lvl         Gtl      
 5 60         RL       50-100      10000-… None  Lvl         Gtl      
 6 50         RL       50-100      10000-… None  Lvl         Gtl      
 7 20         RL       50-100      10000-… None  Lvl         Gtl      
 8 60         RL       50-100      10000-… None  Lvl         Gtl      
 9 50         RM       50-100      <10000  None  Lvl         Gtl      
10 190        RL       <50         <10000  None  Lvl         Gtl      
# … with 1,450 more rows, and 32 more variables: Neighborhood <fct>,
#   BldgType <fct>, HouseStyle <fct>, OverallCond <fct>,
#   YearBuilt <fct>, RoofStyle <fct>, Exterior1st <fct>,
#   ExterCond <fct>, BsmtQual <fct>, BsmtCond <fct>, HeatingQC <fct>,
#   CentralAir <fct>, Electrical <fct>, `1stFlrSF` <fct>,
#   `2ndFlrSF` <fct>, GrLivArea <fct>, KitchenQual <fct>,
#   Functional <fct>, FireplaceQu <fct>, GarageCars <fct>, …

Fase 4.1: Modelo y Flujo (KNN en modo regresión)

Una vez definida la receta, definimos el modelo y lo unimos con nuestra receta creando un flujo de clasificación. Se trata de un nuevo modelo con tune() para poder definir posteriormente los parámetros de forma manual. Activamos, además, el modo regresión con mode = "regression"

# Modelo tuneado
knn_model_tune <-
  nearest_neighbor(mode = "regression", neighbors = tune("k"),
                   weight_func = tune("weight"), dist_power = tune("dist")) |>
  set_engine("kknn")

# Flujo de trabajo
house_wflow_knn <-
  workflow() |>
  add_recipe(house_rec_knn) |>
  add_model(knn_model_tune)

Fase 5.1: Evaluación en Validación Cruzada V-Folds (KNN en modo regresión)

Para esta práctica evaluaremos directamente a través del método de Validación Cruzada V-Folds, que nos proporcionará además una métrica media con la desviación típica para cada modelo.

Creación del grid expandido

Para ello, en primer lugar crearemos un grid con la función expand_grid que nos permitirá definir manualmente los parámetros que queramos sin estar condicionados por una función automática. De esta manera, podremos probar todas las combinaciones que queramos. Tras comprobar varias combinaciones, finalmente nos decantamos por la siguiente.

grid_knn <-
  expand_grid("k" = c(8, 20, 60, 80, 100),
              "weight" = c("inv", "gaussian"),
              "dist" = c(0.01, 0.9, 2, 5, 15, 20))
grid_knn
# A tibble: 60 × 3
       k weight    dist
   <dbl> <chr>    <dbl>
 1     8 inv       0.01
 2     8 inv       0.9 
 3     8 inv       2   
 4     8 inv       5   
 5     8 inv      15   
 6     8 inv      20   
 7     8 gaussian  0.01
 8     8 gaussian  0.9 
 9     8 gaussian  2   
10     8 gaussian  5   
# … with 50 more rows

Aplicación del proceso de validación

Aplicamos en esta fase el proceso de Validación Cruzada V-Folds.

house_cv_folds <- vfold_cv(data = house_train, v = 8, repeats = 4)
house_cv_folds
#  8-fold cross-validation repeated 4 times 
# A tibble: 32 × 3
   splits             id      id2  
   <list>             <chr>   <chr>
 1 <split [1277/183]> Repeat1 Fold1
 2 <split [1277/183]> Repeat1 Fold2
 3 <split [1277/183]> Repeat1 Fold3
 4 <split [1277/183]> Repeat1 Fold4
 5 <split [1278/182]> Repeat1 Fold5
 6 <split [1278/182]> Repeat1 Fold6
 7 <split [1278/182]> Repeat1 Fold7
 8 <split [1278/182]> Repeat1 Fold8
 9 <split [1277/183]> Repeat2 Fold1
10 <split [1277/183]> Repeat2 Fold2
# … with 22 more rows
set.seed(05492)

clusters <- detectCores() - 1
make_cluster <- makeCluster(clusters)
registerDoParallel(make_cluster)

house_knn_vf <-
 house_wflow_knn |> 
  tune_grid(resamples = house_cv_folds,
            grid = grid_knn,
            control = control_grid(verbose = TRUE, allow_par = TRUE,
                                   pkgs = c("outliers", "tidyverse")))

# Finalizamos clusters
stopCluster(make_cluster)
registerDoSEQ()

house_knn_vf |> collect_metrics()
# A tibble: 120 × 9
       k weight  dist .metric .estimator    mean     n std_err .config
   <dbl> <chr>  <dbl> <chr>   <chr>        <dbl> <int>   <dbl> <chr>  
 1     8 gauss…  0.01 rmse    standard   3.83e+4    32 1.19e+3 Prepro…
 2     8 gauss…  0.01 rsq     standard   7.71e-1    32 9.62e-3 Prepro…
 3    20 gauss…  0.01 rmse    standard   3.93e+4    32 1.29e+3 Prepro…
 4    20 gauss…  0.01 rsq     standard   7.64e-1    32 1.03e-2 Prepro…
 5    60 gauss…  0.01 rmse    standard   4.23e+4    32 1.36e+3 Prepro…
 6    60 gauss…  0.01 rsq     standard   7.39e-1    32 1.03e-2 Prepro…
 7    80 gauss…  0.01 rmse    standard   4.34e+4    32 1.37e+3 Prepro…
 8    80 gauss…  0.01 rsq     standard   7.30e-1    32 9.94e-3 Prepro…
 9   100 gauss…  0.01 rmse    standard   4.45e+4    32 1.39e+3 Prepro…
10   100 gauss…  0.01 rsq     standard   7.19e-1    32 9.94e-3 Prepro…
# … with 110 more rows

Elección del mejor de los modelos

Con show_best() se muestran a continuación los mejores modelos según `rsq``.

# Se muestran los mejores según rsq
house_knn_vf |> show_best("rsq")
# A tibble: 5 × 9
      k weight    dist .metric .estimator  mean     n std_err .config 
  <dbl> <chr>    <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>   
1     8 gaussian   0.9 rsq     standard   0.823    32 0.00774 Preproc…
2     8 inv        0.9 rsq     standard   0.820    32 0.00802 Preproc…
3    20 gaussian   0.9 rsq     standard   0.817    32 0.00825 Preproc…
4    20 inv        0.9 rsq     standard   0.810    32 0.00872 Preproc…
5    60 gaussian   0.9 rsq     standard   0.803    32 0.00825 Preproc…

Finalización del flujo

Tras elegir el mejor, finalizamos el flujo con ese modelo elegido empleando la función finalize_workflow().

# Finalizamos flujo con el mejor modelo según rsq
best_house_knn_vf <- 
  house_knn_vf |> select_best("rsq")

final_wflow_knn <- 
  house_wflow_knn |> 
  finalize_workflow(best_house_knn_vf)

final_wflow_knn
══ Workflow ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: nearest_neighbor()

── Preprocessor ──────────────────────────────────────────────────────
28 Recipe Steps

• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• ...
• and 18 more steps.

── Model ─────────────────────────────────────────────────────────────
K-Nearest Neighbor Model Specification (regression)

Main Arguments:
  neighbors = 8
  weight_func = gaussian
  dist_power = 0.9

Computational engine: kknn

Fase 4.2: Modelo y Flujo (Árboles en modo regresión) con poda

A continuación, procederemos con el modelo para la receta de árboles. El modelo consistirá en un árbol en modo regresión que implementaremos con decision_tree() y mode = "regression". De nuevo, crearemos el modelo con tune() para poder definir posteriormente los parámetros de forma manual. Generamos un flujo de trabajo en función del criterios de impureza de Gini.

Por otro lado, añadiremos también directamente el parámetro cost_complexity, que se usa en la fase de poda (prune), y que introduce un componente penalizador para evitar modelos sobreajustados.

# Modelo con tres parámetros tuneados
decision_tree_gini <-
  decision_tree(tree_depth = tune("depth"), min_n = tune("min_n"), 
                cost_complexity = tune("cost")) |> 
  set_engine("rpart") |> 
  set_mode("regression")
# Flujo de trabajo
house_tree_gini_wflow <-
  workflow() |>
  add_recipe(house_rec_arbol) |>
  add_model(decision_tree_gini)

Fase 5.2: Evaluación en Validación Cruzada V-Folds (Árboles en modo regresión)

Para este modelo, volvemos a evaluar a través del método de Validación Cruzada V-Folds, que nos proporcionará además una métrica media con la desviación típica para cada modelo.

Creación del grid expandido

Para ello, en primer lugar crearemos un grid con la función expand_grid que nos permitirán definir manualmente los parámetros que queramos sin estar condicionados por una función automática. De nuevo, tras comprobar varias combinaciones, nos decantaremos finalmente por el siguiente.

grid_tree_gini <-
  expand_grid("depth" = c(4, 6, 7, 8, 10),
              "min_n" = c(45, 50, 55, 100, 150),
              "cost" = 10^c(-5, -4, -3, -0.5))

Crearemos en este caso 100 modelos (\(5*5*4\)).

Aplicación del proceso de validación

Aplicamos en esta fase el proceso de Validación Cruzada V-Folds.

set.seed(05492)

clusters <- detectCores() - 1
make_cluster <- makeCluster(clusters)
registerDoParallel(make_cluster)

house_tree_gini <-
  house_tree_gini_wflow |> 
  tune_grid(resamples = house_cv_folds,
            grid = grid_tree_gini,
            control = control_grid(verbose = TRUE, allow_par = TRUE,
                                   pkgs = c("outliers", "tidyverse"), save_pred = TRUE))

# Finalizamos clusters
stopCluster(make_cluster)
registerDoSEQ()

house_tree_gini |> collect_metrics()
# A tibble: 200 × 9
      cost depth min_n .metric .estimator      mean     n   std_err
     <dbl> <dbl> <dbl> <chr>   <chr>          <dbl> <int>     <dbl>
 1 0.00001     4    45 rmse    standard   47332.       32 1049.    
 2 0.00001     4    45 rsq     standard       0.645    32    0.0110
 3 0.0001      4    45 rmse    standard   47332.       32 1049.    
 4 0.0001      4    45 rsq     standard       0.645    32    0.0110
 5 0.001       4    45 rmse    standard   47331.       32 1049.    
 6 0.001       4    45 rsq     standard       0.645    32    0.0110
 7 0.316       4    45 rmse    standard   62793.       32 1145.    
 8 0.316       4    45 rsq     standard       0.372    32    0.0140
 9 0.00001     4    50 rmse    standard   46869.       32 1058.    
10 0.00001     4    50 rsq     standard       0.651    32    0.0107
# … with 190 more rows, and 1 more variable: .config <chr>

Elección del mejor de los modelos

Como sucedía en el modelo KNN, procedemos a ver los mejores con show_best().

# Se muestran los mejores según rsq
house_tree_gini |>  show_best("rsq", n = 10)
# A tibble: 10 × 9
      cost depth min_n .metric .estimator  mean     n std_err .config 
     <dbl> <dbl> <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>   
 1 0.00001     8    50 rsq     standard   0.699    32  0.0105 Preproc…
 2 0.0001      8    50 rsq     standard   0.699    32  0.0105 Preproc…
 3 0.001       8    50 rsq     standard   0.699    32  0.0106 Preproc…
 4 0.001      10    50 rsq     standard   0.698    32  0.0107 Preproc…
 5 0.00001    10    50 rsq     standard   0.698    32  0.0107 Preproc…
 6 0.0001     10    50 rsq     standard   0.698    32  0.0107 Preproc…
 7 0.00001     8    55 rsq     standard   0.697    32  0.0104 Preproc…
 8 0.0001      8    55 rsq     standard   0.697    32  0.0104 Preproc…
 9 0.0001     10    55 rsq     standard   0.696    32  0.0106 Preproc…
10 0.00001    10    55 rsq     standard   0.696    32  0.0106 Preproc…

Finalización del flujo

Tras elegir el mejor, finalizamos el flujo con ese modelo elegido empleando la función finalize_workflow().

# Finalizamos flujo con el mejor modelo según rsq
best_house_tree_gini <-
  house_tree_gini |> select_best("rsq")

final_wf_tree_gini <- 
  house_tree_gini_wflow |> 
  finalize_workflow(best_house_tree_gini)

final_wf_tree_gini
══ Workflow ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ──────────────────────────────────────────────────────
24 Recipe Steps

• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• ...
• and 14 more steps.

── Model ─────────────────────────────────────────────────────────────
Decision Tree Model Specification (regression)

Main Arguments:
  cost_complexity = 1e-05
  tree_depth = 8
  min_n = 50

Computational engine: rpart

Comparativa final entre modelos: KNN y árbol en modo regresión

A continuación un resumen de sus métricas para hacer la comparativa:

Métricas de los modelos

# Modelo KNN
house_knn_vf |> 
  collect_metrics() |> 
  filter(.config == pull(select(best_house_knn_vf, .config)))
# A tibble: 2 × 9
      k weight   dist .metric .estimator    mean     n std_err .config
  <dbl> <chr>   <dbl> <chr>   <chr>        <dbl> <int>   <dbl> <chr>  
1     8 gaussi…   0.9 rmse    standard   3.40e+4    32 1.06e+3 Prepro…
2     8 gaussi…   0.9 rsq     standard   8.23e-1    32 7.74e-3 Prepro…
# Modelo Árbol Gini
house_tree_gini |> 
  collect_metrics() |> 
  filter(.config == pull(select(best_house_tree_gini, .config)))
# A tibble: 2 × 9
     cost depth min_n .metric .estimator    mean     n std_err .config
    <dbl> <dbl> <dbl> <chr>   <chr>        <dbl> <int>   <dbl> <chr>  
1 0.00001     8    50 rmse    standard   4.36e+4    32 1.16e+3 Prepro…
2 0.00001     8    50 rsq     standard   6.99e-1    32 1.05e-2 Prepro…

Tanto en RMSE como en RSQ el modelo KNN parece mejor: presenta un menor RMSE y también un mayor RSQ.

Predicción para el mejor modelo: KNN

Nos quedamos finalmente con el KNN.

# Ajuste
house_knn_fit <- 
  fit(final_wflow_knn, house_train)

# Predicciones
pred_knn <-
  predict(house_knn_fit, house_test)
summary(pred_knn)

# Visualización de las predicciones para cada `Id`
final_pred <- 
  data.frame(Id = c(1461:2919), pred_knn) |> 
  dplyr::rename(SalePrice = .pred)

head(final_pred)

# Exportamos nuestro dataframe para subirlo a Kaggle
write.csv(final_pred, file = 'house_entrega_knn.csv', row.names = FALSE)

Fase 4.3: Receta, Modelo y Flujo (Regresión univariante)

El tercer modelo que vamos a probar para predecir será la regresión univariante. En este tipo de regresión, seleccionaremos la variable que mayor correlación presenta respecto de nuestra variable objetivo SalePrice. A continuación, visualizaremos de nuevo la matriz de correlaciones y seleccionaremos la susodicha variable.

Preliminares

Selección de nuestra variable (matriz de correlaciones)

cor_matrix <- 
  house_train |> select(where(is.numeric)) |> cor(use = "pairwise.complete.obs", method = "pearson")
cor_matrix |>
  corrplot(method = "number", tl.cex = 0.55, number.cex = 0.7, type = "lower")

A la vista de la matriz, la variable que mayor correlación presenta respecto de nuestra objetivo es GrLivArea con un coeficiente de correlación del 0.71. De esta manera, pasaremos a predecir nuestra objetivo SalePrice en función de GrLivArea.

Comprobación de datos ausentes

Como la regresión univariante no admite ausentes ni en la predictora ni en la objetivo, lo que vamos a hacer es utilizar únicamente el conjunto de train en el que nuestra variable objetivo está completamente etiquetada.

ausentes <- 
  apply(house_train, 2, function(x) sum(is.na(x)))

ausentes_tb <- 
  tibble(Variable = names(house_train), Ausentes = ausentes) |> 
  filter(Variable == "GrLivArea" | Variable == "SalePrice")
ausentes_tb
# A tibble: 2 × 2
  Variable  Ausentes
  <chr>        <int>
1 GrLivArea        0
2 SalePrice        0

Como se puede observar, no hay ausentes ni en nuestra predictora ni en nuestra objetivo.

Comprobación de outliers

house[!is.na(house$SalePrice), ] |> 
  ggplot(aes(x = GrLivArea, y = SalePrice)) +
  geom_point(col = "#EB9891") + 
  geom_smooth(method = "lm", se = FALSE, color = "black", aes(group = 1)) +
  scale_y_log10(breaks= seq(0, 800000, by = 100000), 
                     labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_x_log10(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_text_repel(aes(label = ifelse(house$GrLivArea[!is.na(house$SalePrice)] > 4000, 
                                     rownames(house), ""))) +
  geom_text_repel(aes(label = ifelse(house$SalePrice[!is.na(house$SalePrice)] < 50000, 
                                     rownames(house), ""))) +
  theme_minimal()

Si se grafica la relación siguiendo una escala logarítmica, podremos visualizar claramente esa correlación directa entre variables: a mayor superficie habitable, mayor precio. Además, podemos visualizar también varios valores atípicos. Respecto de la variableGrLivArea, tenemos las tuplas correspondientes a los Id 692, 1183, 524 y 1299. Pasaremos ahora a detectarlos numéricamente, pero antes transformaremos nuestras variables a logarítmicas.

house_log <-
  house_train |> 
  mutate(SalePrice = log(SalePrice),
         GrLivArea = log(GrLivArea))
cor_matrix <- 
  house_log |> 
  select(where(is.numeric)) |>  cor(use = "pairwise.complete.obs", method = "pearson")

cor_matrix |>
  corrplot(method = "number", tl.cex = 0.55, number.cex = 0.7, type = "lower")

El coeficiente de correlación ha aumentado ligeramente, de 0.71 a 0.73.

gg1 <- ggplot(house_train, aes(sample = SalePrice))+ 
  stat_qq()+
  stat_qq_line()+ 
  labs(x = "SalePrice (lineal)", y = "") +
  theme_minimal()

gg2 <- ggplot(house_train, aes(sample = GrLivArea))+ 
  stat_qq()+
  stat_qq_line()+ 
  labs(x = "GrLivArea (lineal)", y = "") +
  theme_minimal()

gg3 <- ggplot(house_log, aes(sample = SalePrice))+ 
  stat_qq()+
  stat_qq_line()+ 
  labs(x = "SalePrice (logarítmica)", y = "") +
  theme_minimal()

gg4 <- ggplot(house_log, aes(sample = GrLivArea))+ 
  stat_qq()+
  stat_qq_line()+ 
  labs(x = "GrLivArea (logarítmica)", y = "") +
  theme_minimal()

multiplot(gg1, gg2, gg3, gg4, cols = 2)

En ambas variables podemos observar como la transformación a logaritmo favorece la relación lineal. Ahora sí, pasemos a detectar outliers sobre nuestra predictora:

house_log |> 
  ggplot(aes(x = SalePrice)) +
  geom_density(alpha = .8, fill="#EB9891") +
  labs(title = "Distribución del precio de las viviendas (escala logarítmica)", x = "Precio", y = NULL) +
  theme_minimal() +
  theme(text = element_text(face = "bold", size = 15), plot.title = element_text(hjust = 0.5), 
        axis.title.x = element_text(vjust=-0.5)) + 
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_y_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_vline(aes(xintercept = mean(SalePrice, na.rm = T), linetype = "Media"), colour = "black", size = .8) +
  geom_vline(aes(xintercept = median(SalePrice, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) + 
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted"))
house_log |> 
  ggplot(aes(x = GrLivArea)) +
  geom_density(alpha = .8, fill="#EB9891") +
  labs(title = "Distribución de la superficie habitable de las viviendas en pies cuadrados (escala logarítmica)", x = "Pies cuadrados", y = NULL) +
  theme_minimal() +
  theme(text = element_text(face = "bold", size = 15), plot.title = element_text(hjust = 0.5), 
        axis.title.x = element_text(vjust=-0.5)) + 
  scale_x_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  scale_y_continuous(labels = comma_format(big.mark = " ", decimal.mark = ",")) +
  geom_vline(aes(xintercept = mean(GrLivArea, na.rm = T), linetype = "Media"), colour = "black", size = .8) +
  geom_vline(aes(xintercept = median(GrLivArea, na.rm = T), linetype = "Mediana"), colour = "black", size = .8) + 
  scale_linetype_manual(name = "Medidas", values = c(Media = "solid", Mediana = "dotted"))

Ambas variables parecen seguir distribuciones más o menos simétricas, por lo que detectaremos los outliers de la predictora por la media.

house_log |>  
  mutate(outliers = abs(scores(GrLivArea, type = "z")) > 2.5) |> 
  filter(outliers)
# A tibble: 15 × 40
   MSSubClass MSZoning LotFrontage LotArea Alley LandContour LandSlope
   <fct>      <fct>          <dbl>   <dbl> <fct> <fct>       <fct>    
 1 30         RM                60    6324 None  Lvl         Gtl      
 2 75         RM                90   22950 None  Lvl         Gtl      
 3 75         RM                87   18386 None  Lvl         Gtl      
 4 60         RL               130   40094 None  Bnk         Gtl      
 5 30         RL                58    9098 None  Lvl         Gtl      
 6 20         RL                50    5000 None  Low         Mod      
 7 190        RH                60   10896 Pave  Bnk         Gtl      
 8 60         RL               104   21535 None  Lvl         Gtl      
 9 30         RM                50    6000 None  Lvl         Gtl      
10 20         C (all)           50    9000 None  Lvl         Gtl      
11 30         RL                60    8400 None  Bnk         Gtl      
12 60         RL               118   35760 None  Lvl         Gtl      
13 60         RL               160   15623 None  Lvl         Gtl      
14 50         RL                NA   14100 None  Lvl         Mod      
15 60         RL               313   63887 None  Bnk         Gtl      
# … with 33 more variables: Neighborhood <fct>, BldgType <fct>,
#   HouseStyle <fct>, OverallCond <dbl>, YearBuilt <dbl>,
#   RoofStyle <fct>, Exterior1st <fct>, ExterCond <int>,
#   BsmtQual <int>, BsmtCond <int>, HeatingQC <int>,
#   CentralAir <fct>, Electrical <fct>, `1stFlrSF` <dbl>,
#   `2ndFlrSF` <dbl>, GrLivArea <dbl>, KitchenQual <dbl>,
#   Functional <fct>, FireplaceQu <int>, GarageCars <dbl>, …

Durante la creación de la receta le imputaremos a estos outliers la media.

Receta para el método de regresión univariante

house_rec_uni <-
  recipe(data = house_train, SalePrice ~ GrLivArea) |> 
  step_log(SalePrice, base = 10) |> 
  step_log(GrLivArea, base = 10) |> 
  step_mutate(GrLivArea = ifelse(abs(scores(GrLivArea, type = "z")) > 2.5, NA, GrLivArea)) |> 
  step_impute_mean(GrLivArea)
house_rec_uni
Recipe

Inputs:

      role #variables
   outcome          1
 predictor          1

Operations:

Log transformation on SalePrice
Log transformation on GrLivArea
Variable mutation for ifelse(abs(scores(GrLivArea, typ...
Mean imputation for GrLivArea

Horneado

Por último, horneamos nuestra receta para comprobar que todo funciona correctamente.

bake(house_rec_uni |> prep(), new_data = NULL) 
# A tibble: 1,460 × 2
   GrLivArea SalePrice
       <dbl>     <dbl>
 1      3.23      5.32
 2      3.10      5.26
 3      3.25      5.35
 4      3.23      5.15
 5      3.34      5.40
 6      3.13      5.16
 7      3.23      5.49
 8      3.32      5.30
 9      3.25      5.11
10      3.03      5.07
# … with 1,450 more rows

Modelo y Flujo

# Creación del modelo
reg_lineal <- 
  linear_reg() |> set_mode("regression") |> set_engine("lm")
# Creación del flujo
reg_wflow_uni <-
  workflow() |> 
  add_model(reg_lineal) |> 
  add_recipe(house_rec_uni)
reg_wflow_uni
══ Workflow ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ──────────────────────────────────────────────────────
4 Recipe Steps

• step_log()
• step_log()
• step_mutate()
• step_impute_mean()

── Model ─────────────────────────────────────────────────────────────
Linear Regression Model Specification (regression)

Computational engine: lm

Ajuste

set.seed(05492)

reg_fit_uni <-
  reg_wflow_uni |> fit(data = house_train)
reg_fit_uni 
══ Workflow [trained] ════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ──────────────────────────────────────────────────────
4 Recipe Steps

• step_log()
• step_log()
• step_mutate()
• step_impute_mean()

── Model ─────────────────────────────────────────────────────────────

Call:
stats::lm(formula = ..y ~ ., data = data)

Coefficients:
(Intercept)    GrLivArea  
     2.4439       0.8804  
# Resumen del ajuste
tidy(reg_fit_uni)
# A tibble: 2 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    2.44     0.0752      32.5 1.27e-174
2 GrLivArea      0.880    0.0238      37.0 1.42e-211

Como se puede observar en las métricas, el modelo que se ha decidido crear se trata de un modelo doble logarítmico. La interpretación de \(\beta_1\) es la siguiente: si el logaritmo de la superficie habitable de una vivienda aumenta en un 1 % se estima que, en media, el logaritmo del precio de la vivienda se incrementará en un 0,88 %. Esta relación no es más que la elasticidad entre la superficie habitable de una vivienda y su precio y se trata, además, de una relación constante. Por otro lado, \(\beta_0\) puede interpretarse como el valor medio del precio de la vivienda cuando la superficie habitable toma valor cero. En este caso, \(\beta_0\) toma el valor de 2.44, aunque en muchas situaciones su interpretación carece de sentido.

Por otro lado, en la salida del modelo podemos observar como el p-valor coligado al parámetro \(\beta_1\) y al estadístico para este contraste es p = 1.415987e-211. Si trabajamos con un nivel de significación habitual del 5 % (0.05), el p-valor es menor que ese nivel de significación, por lo que podemos rechazar la hipótesis nula, esto es, la no significación individual de la variable log(GrLivArea) al 5 % de significación. Por ende, a la luz de los resultados, \(\beta_1\) es significativo individualmente, o lo que es lo mismo: el efecto de la variable log(GrLivArea) depende de la variable log(SalePrice) al 5 % de nivel de significación.

Intervalo de confianza e índices performáticos del modelo

# Intervalo de confianza al 95 %
confint(reg_fit_uni |> extract_fit_engine())
                2.5 %   97.5 %
(Intercept) 2.2963505 2.591491
GrLivArea   0.8336623 0.927106

Los intervalos de confianza nos ofrecen un rango de valores para \(\beta_j\) en donde no se rechazaría la hipótesis \(H_0\). En este caso concreto, el intervalo de confianza nos proporcionará un rango de valores para la variable explicativa log(GrLivArea) en donde podríamos comprobar su significación individual dentro del modelo.

reg_fit_uni |> extract_fit_engine() |>  performance()
# Indices of model performance

AIC       |      AICc |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------------------
-1931.953 | -1931.937 | -1916.095 | 0.484 |     0.483 | 0.125 | 0.125

Diagnosis

check_model(reg_fit_uni |> extract_fit_engine())

Supuesto de linealidad

En el gráfico, la línea es bastante atendencial, se ve bastante monótona en todo su recorrido. Comprobémoslo analíticamente con un ANOVA entre residuos y predictora:

adjustment <- 
  reg_fit_uni |>  extract_fit_engine()
lm(adjustment$residuals ~ adjustment$fitted.values) |>  anova()
Analysis of Variance Table

Response: adjustment$residuals
                           Df Sum Sq  Mean Sq F value Pr(>F)
adjustment$fitted.values    1  0.000 0.000000       0      1
Residuals                1458 22.668 0.015547               
lm(adjustment$residuals ~ I(adjustment$fitted.values^2) + adjustment$fitted.values) |>  anova()
Analysis of Variance Table

Response: adjustment$residuals
                                Df  Sum Sq  Mean Sq F value Pr(>F)
I(adjustment$fitted.values^2)    1  0.0000 0.000007  0.0004 0.9831
adjustment$fitted.values         1  0.0332 0.033153  2.1340 0.1443
Residuals                     1457 22.6347 0.015535

Como se puede observar, el p-valor de la prueba F de Fisher-Snedecor es igual a 1, por lo que, a un nivel de significación del 0.05, no podemos rechazar la hipótesis nula, esto es, la no existencia de una relación lineal entre la predictora (GrLivArea) y sus residuos. Por tanto, ello implica que no parece existir tendencia lineal entre residuos y predictora en nuestro modelo. Lo mismo ocurre en el segundo caso: tampoco parece existir tendencia cuadrática entre residuos y predictora.

Supuesto de homocedasticidad

En el gráfico, la línea es bastante atendencial, se ve bastante monótona en todo su recorrido (excepto quizá en los extremos). Comprobémoslo analíticamente con un test de heterocedasticidad:

check_heteroscedasticity(adjustment)
Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

Según el test de heterocedasticidad, el supuesto de homocedasticidad no se cumple para nuestro modelo. Esto, muy a menudo, resulta poco realista, pues el supuesto implica que la variabilidad del error de la variable objetivo es la misma para cualquier nivel de nuestra predictora.

ggplot(
  tibble("Observaciones" = 1:length(adjustment$residuals),
         "Residuos" = adjustment$residuals),
  aes(x = Observaciones, y = Residuos)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

En el gráfico, la recta de regresión se aprecia constante, en torno a 0, y los residuos se agrupan en torno a una banda también relativamente constante.

Supuesto de normalidad de los residuos

En el gráfico, los percentiles empíricos de los residuos se acoplan bastante bien a la diagonal, excepto en los extremos. Comprobémoslo analíticamente con un contraste de normalidad:

library(olsrr)
ols_test_normality(adjustment$residuals)
-----------------------------------------------
       Test             Statistic       pvalue  
-----------------------------------------------
Shapiro-Wilk              0.9642         0.0000 
Kolmogorov-Smirnov        0.0685         0.0000 
Cramer-von Mises         383.4014        0.0000 
Anderson-Darling         11.4357         0.0000 
-----------------------------------------------

Como se puede observar, los p-valor de los distintos test son iguales a 0, por lo que, a un nivel de significación del 0.05, podemos rechazar la hipótesis de normalidad. Por tanto, no podemos asumir a normalidad en nuestro modelo. Se han hecho pruebas con cambios en la receta añadiendo distintas transformaciones a nuestra objetivo (step_YeoJohnson(), step_BoxCox(), step_sqrt()), pero tampoco ha sido posible.

Incorrelación de los residuos

library(car)
durbinWatsonTest(adjustment)
 lag Autocorrelation D-W Statistic p-value
   1      0.01831851      1.963326   0.488
 Alternative hypothesis: rho != 0

Como se puede observar, el p-valor de la prueba Durbin-Watson es igual a 0.476, por lo que, a un nivel de significación del 0.05, no podemos rechazar la hipótesis nula, esto es, la incorrelación de los residuos.

Influencia de los residuos

En el gráfico original de la diagnosis, todas las observaciones parecen estar contenidas entre las líneas punteadas al 0.5.

Fase 5.3: Evaluación del modelo (Regresión univariante)

Resumen

glance(reg_fit_uni)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic   p.value    df logLik
      <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl>
1     0.484         0.483 0.125     1366. 1.42e-211     1   969.
# … with 5 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
#   df.residual <int>, nobs <int>

Nuestro \(R^2\) es igual a 0.4838, por lo que el ratio de información explicada del modelo es del 48.38 %. Al ser un modelo tan simple con una única variable, la interpretación del resto de parámetros no tiene mucho sentido. Los dejaremos para los próximos modelos. A continuación, predeciremos y evaluaremos en un nuevo conjunto test que extraeremos de house_train, puesto que nuestro conjunto house_test original tiene la variable objetivo SalePrice repleta de NA.

Evaluación y predicciones en test

# Hacemos el split. Lo hacemos del 0.6 ya que hay pocos datos en train. 
set.seed(05492)

split_house <- 
  initial_split(house_train, prop = 0.6)

# Predecimos en test
fit_reg_uni <- 
  reg_wflow_uni |> last_fit(split = split_house)

# Evaluamos en test
fit_reg_uni |>  collect_metrics()
# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       0.126 Preprocessor1_Model1
2 rsq     standard       0.455 Preprocessor1_Model1

Visualización de errores en test

# Errores en test
pred_test <-
  fit_reg_uni |>
  collect_predictions() |>
  mutate(error = SalePrice - .pred)
pred_test
# A tibble: 584 × 6
   id               .pred  .row SalePrice .config                error
   <chr>            <dbl> <int>     <dbl> <chr>                  <dbl>
 1 train/test split  5.17     2      5.26 Preprocessor1_Model1  0.0874
 2 train/test split  5.31     3      5.35 Preprocessor1_Model1  0.0394
 3 train/test split  5.29     7      5.49 Preprocessor1_Model1  0.198 
 4 train/test split  5.37     8      5.30 Preprocessor1_Model1 -0.0715
 5 train/test split  5.31     9      5.11 Preprocessor1_Model1 -0.194 
 6 train/test split  5.09    11      5.11 Preprocessor1_Model1  0.0179
 7 train/test split  5.04    13      5.16 Preprocessor1_Model1  0.116 
 8 train/test split  5.17    15      5.20 Preprocessor1_Model1  0.0273
 9 train/test split  5.08    17      5.17 Preprocessor1_Model1  0.0929
10 train/test split  5.20    20      5.14 Preprocessor1_Model1 -0.0520
# … with 574 more rows

En esta tabla se puede observar cómo, en nuestro modelo doble logarítmico, las estimaciones se ajustan en gran medida a los valores reales de la objetivo SalePrice. Visualizaremos ahora estos valores a través de un gráfico.

g1 <- pred_test |> 
  ggplot(mapping = aes(x = .pred, y = SalePrice)) +
  geom_point(color = "#56BCC2", alpha = 0.6, size = 4) +
  geom_abline(intercept = 0, slope = 1, color = "#EB9891", size = 1.2) +
  theme_minimal() + 
  labs(title = "Resultados de la regresión lineal univariante",
       subtitle = "Los valores predichos deberían estar cercanos a la diagonal",
       x = "Predicciones",
       y = "Valores reales")

g2 <- pred_test |> 
  select(.pred, SalePrice) |> 
  gather(Distribución, value) |> 
  ggplot(aes(x = value, color = Distribución, fill = Distribución)) + 
  geom_density(alpha = 0.6) + 
  theme_minimal() + 
  labs(title = "Distribución de las predicciones sobre los valores reales de SalePrice",
       x = "Distribuciones",
       y = "Frecuencia")

multiplot(g1, g2)

Como se puede apreciar en los gráficos, las predicciones son un poco reguleras, principalmente por el reducido tamaño muestral y por el bajo \(R^2\) (finalmente, el ratio de información explicada del modelo fue del 49.88 %). Tampoco se han cumplido varias de las hipótesis. Desde una regresión lineal univariante, en la que únicamente disponemos de una única variable para predecir la objetivo, poco más se puede hacer. En los siguientes apartados avanzaremos hacia la regresión multivariante con la creación de dos nuevos modelos.

Fase 4.4: Receta, Modelo y Flujo (Regresión multivariante: modelo saturado)

Para este nuevo modelo de regresión vamos a emplear todas las variables. Volvemos a generar la receta, pero esta vez enfrentando nuestra objetivo al total de variables predictoras que tenemos en el dataset.

Receta para el método de regresión multivariante

Aplicación de roles

Volvemos a definir una nueva receta indicándole el conjunto donde tenemos validación y train, y enfrentando nuestra variable objetivo SalePrice a todas las demás. Después, asignamos posibles roles, sujetos a modificación, que nos permitan diferenciar acciones entre las variables (sobre todo en la sección outliers).

# Receta
house_rec_multi <-
  # Fórmula y datos
  recipe(data = house_train, SalePrice ~ .)|>
  # Roles
  add_role(where(is.factor), 
           new_role = "cualitativa") |> 
  add_role(where(is.numeric), 
           new_role = "cuantitativa") |> 
  add_role(RemodAdd, New, 
           new_role = "binaria") |> 
  add_role(LotFrontage, LotArea, `1stFlrSF`, `2ndFlrSF`, GrLivArea, Age, 
           new_role = "mediana") |> 
  add_role(all_nominal_predictors(), 
           new_role = "moda") |> 
  add_role(ExterCond, BsmtCond, HeatingQC, GarageQual, KitchenQual,
           new_role = "imputar mediana") |> 
  add_role(OverallCond, BsmtQual, FireplaceQu, TotBath, YearBuilt, GarageCars,
           new_role = "imputar media") |> 
  add_role(SalePrice, GrLivArea,`1stFlrSF`, `2ndFlrSF`,
           new_role = "log")

Reagrupación de las variables cualitativas

En este apartado, reagrupamos las variables cualitativas con step_mutate con lo definido en la fase de exploración. Vamos a reducir el número de categorías para que, al dummyficar, no se creen en exceso sin ser realmente necesarias.

house_rec_multi <- 
  house_rec_multi |> 
  step_mutate(Alley = 
                fct_collapse(Alley, 
                             Grvl_Pave = c("Grvl", "Pave"))) |> 
  step_mutate(LandSlope = 
                fct_collapse(LandSlope, 
                             Mod_Sev = c("Mod", "Sev"))) |>  
  step_mutate(LandContour = 
                fct_collapse(LandContour, 
                             Bnk_HLS_Low = c("Bnk", "HLS", "Low"))) |> 
  step_mutate(RoofStyle = 
                fct_collapse(RoofStyle, 
                             Hip = c("Flat", "Shed", "Mansard", "Hip"),
                             Gam_Gable = c("Gambrel", "Gable"))) |> 
  step_mutate(Exterior1st = 
                fct_collapse(Exterior1st, 
                            `<150k` = c("AsphShn", "BrkComm", "BrkFace",
                                        "CemntBd", "HdBoard", "MetalSd",
                                        "Stucco", "Wd Sdng", "WdShing",
                                        "AsbShng"),
                             `>150k` = c("CBlock", "ImStucc", "Plywood",
                                         "Stone", "VinylSd"))) |> 
  step_mutate(PavedDrive = 
                fct_collapse(PavedDrive, 
                             N_P = c("N", "P"))) |>
  step_mutate(Fence = 
                fct_collapse(Fence, 
                             MnPrv = c("GdWo", "MnWw"),
                             None = c("GdPrv"))) |> 
  step_mutate(MSSubClass = 
                fct_collapse(MSSubClass, 
                             `<150k` = c("30", "40", "45", "50",
                                         "80", "85", "120", "160",
                                         "180", "190"),
                             `>150k` = c("20", "60", "70", "75",
                                         "90", "150"))) |>
  step_mutate(MSZoning = 
                fct_collapse(MSZoning, 
                             RM = c("C (all)", "RH"),
                             RL = c("FV"))) |> 
  step_mutate(Neighborhood = 
                fct_collapse(Neighborhood, 
                             `<125k` = c("BrDale", "BrkSide", "Edwards",
                                         "IDOTRR", "MeadowV", "OldTown"),
                             `125k-150k` = c("Blmngtn", "NAmes", "NPkVill",
                                             "SWISU", "Sawyer"),
                             `150k-200k` = c("Blueste", "CollgCr", "Gilbert",
                                             "Mitchel", "NWAmes", "SawyerW"),
                             `>200k` = c("ClearCr", "Crawfor", "NoRidge",
                                         "NridgHt", "Somerst", "StoneBr",
                                         "Timber", "Veenker"))) |> 
  step_mutate(BldgType = 
                fct_collapse(BldgType, 
                             Fam_TwnhsE = c("1Fam", "TwnhsE"),
                             fmCon_Dup_Twnhs = c("2fmCon", "Duplex", "Twnhs"))) |>  
  step_mutate(HouseStyle = 
                fct_collapse(HouseStyle, 
                             Story_Fin_SLvl = c("2Story", "2.5Fin", "SLvl"),
                             fmCon_Dup_Twnhs = c("1.5Fin", "1.5Unf", "1Story", 
                                                 "2.5Unf", "SFoyer"))) |> 
  step_mutate(Electrical = 
                fct_collapse(Electrical, 
                             Other = c("FuseA", "FuseF", "FuseP", "Mix"))) |>
  step_mutate(Functional = 
                fct_collapse(Functional, 
                             Other = c("Maj1", "Maj2", "Min1", "Min2", 
                                       "Mod", "Sev"))) |>
  step_mutate(SaleType = 
                fct_collapse(SaleType, 
                             Con_CWD_New = c("Con", "CWD", "New"),
                             WD = c("COD", "ConLD", "ConLI", "ConLw",
                                       "Oth"))) |> 
  step_mutate(SaleCondition = 
                fct_collapse(SaleCondition, 
                             Normal_Others = c("Abnorml", "AdjLand", "Alloca",
                                               "Family"))) |>
  step_mutate(MoSold = 
                fct_collapse(MoSold, 
                             Jan_Ap_May_Oct = c("1", "4", "5", "10"),
                             Mar_Jun_Jul = c("3", "6", "7"),
                             Feb_Aug_Sep_Nov_Dec = c("2", "8", "9",
                                                     "11", "12"))) |> 
  step_mutate(YrSold = 
                fct_collapse(YrSold, 
                             `2006-2009` = c("2006", "2007", "2008", "2009")))

Recategorización de las variables cuantitativas

En los métodos de regresión, la recategorización de las variables cuantitativas no es estrictamente necesaria. Por tanto, en este receta, del total de variables cuantitativas no vamos a modificar ninguna.

Outliers

multiplot(box1, box2, box3, box4, box5, box6, box7, box8, box9, box10, cols = 2) 
multiplot(box11, box12, box13, box14, box15, box16, box17, cols = 2)

Si observamos estos gráficos de cajas y bigotes, todas nuestras variables cuantitativas continuas son asimétricas, por lo que se detectarán los outliers y se imputarán los ausentes por la mediana. Para el caso de las semi-cuantitativas (KitchenQual, por ejemplo) se imputarán directamente los ausentes por uno de los dos métodos en función de la simetría o asimetría de la variable. Para el resto de variables cualitativas, imputamos directamente por la moda.

house_rec_multi <-
  house_rec_multi |> 
  # Detección de outliers por la mediana y por la media
  step_mutate(across(has_role("mediana"), function(x) { ifelse(abs(scores(x, type = "mad")) > 3 & !is.na(x), NA, x) })) |>
  # Imputación de ausentes por la mediana, la media y la moda
  step_impute_median(has_role("mediana")) |>
  step_impute_median(has_role("imputar_mediana")) |> 
  step_impute_mean(has_role("imputar_media")) |>
  step_impute_mode(has_role("moda"))

YeoJohnson para la normalidad de los residuos

A fin de que nuestro modelo cumpla los test de normalidad de los residuos, vamos a incluir en la receta la función step_YeoJohnson, que se emplea para eliminar la asimetría en los datos numéricos del modelo.

house_rec_multi <-
  house_rec_multi |> 
  step_YeoJohnson(all_numeric_predictors(), -has_role("log"))

Transformación logarítmica de algunas variables

Transformamos algunas variables a logarítmicas para mejorar su relación.

house_rec_multi <-
  house_rec_multi |> 
  step_log(has_role("log"), offset = 1)

Normalizar por rango

Normalizamos nuestras variables por rango para que todas tengan el mismo peso, entre 0 y 1.

house_rec_multi <-
  house_rec_multi |> 
  step_normalize(all_numeric_predictors())

Variables dummy

Como esta receta es para un modelo de regresión multivariante, debemos dummyficar nuestras variables cualitativas. Para ello, tomamos todas las nominales, menos nuestra variable objetivo.

house_rec_multi <-
  house_rec_multi |>
  step_dummy(all_nominal_predictors())

Filtro de cero varianza

house_rec_multi <-
  house_rec_multi |>
  step_zv(all_predictors())

Filtro de correlación

Aplicamos el filtro de correlación a nuestras variables numéricas.

house_rec_multi <-
  house_rec_multi |> 
  step_corr(all_numeric_predictors(), threshold = 0.6)

Horneado

Por último, horneamos nuestra receta para comprobar que todas nuestras nuevas variables recategorizadas se hayan creado correctamente.

bake(house_rec_multi |>  prep(), new_data = NULL)
# A tibble: 1,460 × 35
   LotArea OverallCond ExterCond BsmtQual BsmtCond HeatingQC
     <dbl>       <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
 1 -0.141       -0.477    -0.194    0.576 -0.00502     0.930
 2  0.106        2.01     -0.194    0.576 -0.00502     0.930
 3  0.414       -0.477    -0.194    0.576 -0.00502     0.930
 4  0.0955      -0.477    -0.194   -0.695  3.56       -0.330
 5  0.877       -0.477    -0.194    0.576 -0.00502     0.930
 6  0.857       -0.477    -0.194    0.576 -0.00502     0.930
 7  0.201       -0.477    -0.194    2.13  -0.00502     0.930
 8  0.257        0.440    -0.194    0.576 -0.00502     0.930
 9 -0.761       -0.477    -0.194   -0.695 -0.00502    -0.330
10 -0.391        0.440    -0.194   -0.695 -0.00502     0.930
# … with 1,450 more rows, and 29 more variables: `1stFlrSF` <dbl>,
#   `2ndFlrSF` <dbl>, GrLivArea <dbl>, FireplaceQu <dbl>,
#   GarageCars <dbl>, GarageQual <dbl>, RemodAdd <dbl>, New <dbl>,
#   Age <dbl>, SalePrice <dbl>, MSSubClass_X.150k <dbl>,
#   MSZoning_RL <dbl>, Alley_None <dbl>, LandContour_Lvl <dbl>,
#   LandSlope_Mod_Sev <dbl>, Neighborhood_X150k.200k <dbl>,
#   Neighborhood_X.200k <dbl>, BldgType_fmCon_Dup_Twnhs <dbl>, …

Flujo y ajuste

# Creación del flujo
reg_wflow_multi <-
  workflow() |> 
  add_model(reg_lineal) |> 
  add_recipe(house_rec_multi)
reg_wflow_multi
══ Workflow ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ──────────────────────────────────────────────────────
29 Recipe Steps

• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• ...
• and 19 more steps.

── Model ─────────────────────────────────────────────────────────────
Linear Regression Model Specification (regression)

Computational engine: lm 
set.seed(05492)

reg_fit_multi <-
  reg_wflow_multi |> fit(data = house_train)
reg_fit_multi
══ Workflow [trained] ════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ──────────────────────────────────────────────────────
29 Recipe Steps

• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• step_mutate()
• ...
• and 19 more steps.

── Model ─────────────────────────────────────────────────────────────

Call:
stats::lm(formula = ..y ~ ., data = data)

Coefficients:
               (Intercept)                     LotArea  
                11.6420832                   0.0186542  
               OverallCond                   ExterCond  
                 0.0466561                  -0.0011050  
                  BsmtQual                    BsmtCond  
                 0.0718343                   0.0044886  
                 HeatingQC                  `1stFlrSF`  
                 0.0216881                   0.0319568  
                `2ndFlrSF`                   GrLivArea  
                -0.0230072                   0.1635733  
               FireplaceQu                  GarageCars  
                 0.0260895                   0.0492290  
                GarageQual                    RemodAdd  
                 0.0161250                  -0.0192380  
                       New                         Age  
                 0.0025791                  -0.0418791  
         MSSubClass_X.150k                 MSZoning_RL  
                -0.0138932                   0.0641965  
                Alley_None             LandContour_Lvl  
                 0.0292293                   0.0263570  
         LandSlope_Mod_Sev     Neighborhood_X150k.200k  
                 0.0253077                   0.0383263  
       Neighborhood_X.200k    BldgType_fmCon_Dup_Twnhs  
                 0.1641738                  -0.0502749  
       RoofStyle_Gam_Gable          Exterior1st_X.150k  
                -0.0199646                  -0.0006763  
              CentralAir_Y            Electrical_SBrkr  
                 0.0688261                   0.0015624  
            Functional_Typ                PavedDrive_Y  
                 0.1148388                   0.0462485  
               Fence_MnPrv  MoSold_Feb_Aug_Sep_Nov_Dec  
                -0.0090422                  -0.0073745  
        MoSold_Mar_Jun_Jul                YrSold_X2010  
                -0.0003882                   0.0058141  
      SaleCondition_Normal  
                 0.0460244  
# Resumen del ajuste
tidy(reg_fit_multi)
# A tibble: 35 × 5
   term        estimate std.error statistic  p.value
   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept) 11.6       0.0351    332.    0       
 2 LotArea      0.0187    0.00521     3.58  3.50e- 4
 3 OverallCond  0.0467    0.00469     9.94  1.49e-22
 4 ExterCond   -0.00110   0.00413    -0.268 7.89e- 1
 5 BsmtQual     0.0718    0.00587    12.2   8.35e-33
 6 BsmtCond     0.00449   0.00428     1.05  2.94e- 1
 7 HeatingQC    0.0217    0.00479     4.52  6.59e- 6
 8 `1stFlrSF`   0.0320    0.0140      2.28  2.27e- 2
 9 `2ndFlrSF`  -0.0230    0.0143     -1.61  1.08e- 1
10 GrLivArea    0.164     0.0166      9.87  2.95e-22
# … with 25 more rows
tidy(reg_fit_multi) |> 
  filter(p.value > 0.5)
# A tibble: 6 × 5
  term                estimate std.error statistic p.value
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl>
1 ExterCond          -0.00110    0.00413   -0.268    0.789
2 New                 0.00258    0.00459    0.561    0.575
3 Exterior1st_X.150k -0.000676   0.00915   -0.0739   0.941
4 Electrical_SBrkr    0.00156    0.0151     0.104    0.917
5 MoSold_Mar_Jun_Jul -0.000388   0.00882   -0.0440   0.965
6 YrSold_X2010        0.00581    0.0117     0.498    0.618

Como se puede observar en las métricas, los p-valor coligados a los parámetros \(\beta\) de las variables ExterCond, New, Exterior1st_X.150k, Electrical_SBrkr, MoSold_Mar_Jun_Jul y YrSold_X2010 son superiores a 0.05. Si trabajamos con un nivel de significación habitual del 5 % (0.05), el p-valor es mayor que ese nivel de significación, por lo que no podemos rechazar la hipótesis nula, esto es, la no significación individual de estas variables al 5 % de significación. Por ende, a la luz de los resultados, ExterCond, New, Exterior1st_X.150k, Electrical_SBrkr, MoSold_Mar_Jun_Jul y YrSold_X2010 son variables no significativas individualmente. Optaremos por eliminarlas posteriormente en la selección de modelos porque solo aportan ruido.

Índices performáticos del modelo

reg_fit_multi |> extract_fit_engine() |>  performance()
# Indices of model performance

AIC       |      AICc |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------------------
-1543.682 | -1541.810 | -1353.379 | 0.879 |     0.876 | 0.139 | 0.141

Nuestro \(R^2\) es igual a 0.879, por lo que el ratio de información explicada del modelo es del 88 %. Este modelo ha duplicado el ratio de información explicada respecto del modelo univariante anterior.

Diagnosis

check_model(reg_fit_multi |> extract_fit_engine())

Supuesto de linealidad

En el gráfico, la línea es bastante atendencial, se ve bastante monótona en todo su recorrido. Comprobémoslo analíticamente con un ANOVA entre residuos y predictora:

adjustment <- 
  reg_fit_multi |>  extract_fit_engine()
lm(adjustment$residuals ~ adjustment$fitted.values) |>  anova()
Analysis of Variance Table

Response: adjustment$residuals
                           Df Sum Sq  Mean Sq F value Pr(>F)
adjustment$fitted.values    1  0.000 0.000000       0      1
Residuals                1458 28.267 0.019387               
lm(adjustment$residuals ~ I(adjustment$fitted.values^2) + adjustment$fitted.values) |>  anova()
Analysis of Variance Table

Response: adjustment$residuals
                                Df  Sum Sq  Mean Sq F value    Pr(>F)
I(adjustment$fitted.values^2)    1  0.0001 0.000118  0.0061 0.9376883
adjustment$fitted.values         1  0.2483 0.248296 12.9119 0.0003374
Residuals                     1457 28.0183 0.019230                  
                                 
I(adjustment$fitted.values^2)    
adjustment$fitted.values      ***
Residuals                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como se puede observar, el p-valor de la prueba F de Fisher-Snedecor es igual a 1, por lo que, a un nivel de significación del 0.05, no podemos rechazar la hipótesis nula, esto es, la no existencia de una relación lineal entre las predictoras y sus residuos. Por tanto, ello implica que no parece existir tendencia lineal entre residuos y predictoras en nuestro modelo. Lo mismo ocurre en el segundo caso: tampoco parece existir tendencia cuadrática entre residuos y predictoras.

Supuesto de homocedasticidad

En el gráfico, la línea sigue una fase decreciente al principio y creciente al final, se ve poco monótona en todo su recorrido. Comprobémoslo analíticamente con un test de heterocedasticidad:

check_heteroscedasticity(adjustment)
Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

Según el test de heterocedasticidad, el supuesto de homocedasticidad no se cumple para nuestro modelo. Esto, muy a menudo, resulta poco realista, pues el supuesto implica que la variabilidad del error de la variable objetivo es la misma para cualquier nivel de nuestra predictora.

ggplot(
  tibble("Observaciones" = 1:length(adjustment$residuals),
         "Residuos" = adjustment$residuals),
  aes(x = Observaciones, y = Residuos)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

En el gráfico, la recta de regresión se aprecia constante, en torno a 0, y los residuos se agrupan en torno a una banda también relativamente constante.

Supuesto de colinearidad

En el gráfico original de la diagnosis, prácticamente todas las observaciones están en la franja verde. Tan solo tres están en la roja (es el mejor resultado que obtuvimos).

Supuesto de normalidad de los residuos

En el gráfico, los percentiles empíricos de los residuos se acoplan bastante bien a la diagonal, excepto en los extremos. Comprobémoslo analíticamente con un contraste de normalidad:

library(olsrr)
ols_test_normality(adjustment$residuals)
-----------------------------------------------
       Test             Statistic       pvalue  
-----------------------------------------------
Shapiro-Wilk              0.9313         0.0000 
Kolmogorov-Smirnov        0.0659         0.0000 
Cramer-von Mises         374.2219        0.0000 
Anderson-Darling         11.2063         0.0000 
-----------------------------------------------

Como se puede observar, los p-valor de los distintos test son iguales a 0, por lo que, a un nivel de significación del 0.05, podemos rechazar la hipótesis de normalidad. Por tanto, no podemos asumir a normalidad en nuestro modelo. Se han hecho pruebas con cambios en la receta añadiendo distintas transformaciones a nuestra objetivo (step_YeoJohnson(), step_BoxCox(), step_sqrt()), pero tampoco ha sido posible.

Incorrelación de los residuos

library(car)
durbinWatsonTest(adjustment)
 lag Autocorrelation D-W Statistic p-value
   1     0.007249536      1.985144   0.828
 Alternative hypothesis: rho != 0

Como se puede observar, el p-valor de la prueba Durbin-Watson es igual a 0.772, por lo que, a un nivel de significación del 0.05, no podemos rechazar la hipótesis nula, esto es, la incorrelación de los residuos.

Influencia de los residuos

En el gráfico original de la diagnosis, todas las observaciones parecen estar contenidas entre las líneas punteadas al 0.5.

Fase 5.4: Evaluación del modelo (Regresión multivariante: modelo saturado)

Resumen

glance(reg_fit_multi)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik    AIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>
1     0.879         0.876 0.141      303.       0    34   808. -1544.
# … with 4 more variables: BIC <dbl>, deviance <dbl>,
#   df.residual <int>, nobs <int>

Nuestro \(R^2\) es igual a 0.8786, por lo que el ratio de información explicada del modelo es del 88 %. Este modelo ha duplicado el ratio de información explicada respecto del modelo univariante anterior.

Evaluación y predicciones en test

# Hacemos el split. Lo hacemos del 0.6 ya que hay pocos datos en train. 
set.seed(05492)

split_house <- 
  initial_split(house_train, prop = 0.6)

# Predecimos en test
reg_fit_multi <- 
  reg_wflow_multi |> last_fit(split = split_house)

# Evaluamos en test
reg_fit_multi |>  collect_metrics()
# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       0.152 Preprocessor1_Model1
2 rsq     standard       0.850 Preprocessor1_Model1

Visualización de errores en test

set.seed(05492)

# Errores en test
pred_test <-
  reg_fit_multi |>
  collect_predictions() |>
  mutate(error = SalePrice - .pred)
pred_test
# A tibble: 584 × 6
   id               .pred  .row SalePrice .config                error
   <chr>            <dbl> <int>     <dbl> <chr>                  <dbl>
 1 train/test split  12.3     2      12.1 Preprocessor1_Mode… -0.232  
 2 train/test split  12.2     3      12.3 Preprocessor1_Mode…  0.112  
 3 train/test split  12.6     7      12.6 Preprocessor1_Mode…  0.0661 
 4 train/test split  12.3     8      12.2 Preprocessor1_Mode… -0.132  
 5 train/test split  11.6     9      11.8 Preprocessor1_Mode…  0.137  
 6 train/test split  11.7    11      11.8 Preprocessor1_Mode…  0.0424 
 7 train/test split  11.6    13      11.9 Preprocessor1_Mode…  0.233  
 8 train/test split  11.8    15      12.0 Preprocessor1_Mode…  0.146  
 9 train/test split  11.9    17      11.9 Preprocessor1_Mode… -0.00848
10 train/test split  11.6    20      11.8 Preprocessor1_Mode…  0.195  
# … with 574 more rows

En esta tabla se puede observar cómo, en nuestro modelo, las estimaciones se ajustan en gran medida a los valores reales de la objetivo SalePrice. Visualizaremos ahora estos valores a través de un gráfico.

g1 <- pred_test |> 
  ggplot(mapping = aes(x = .pred, y = SalePrice)) +
  geom_point(color = "#56BCC2", alpha = 0.6, size = 4) +
  geom_abline(intercept = 0, slope = 1, color = "#EB9891", size = 1.2) +
  theme_minimal() + 
  labs(title = "Resultados de la regresión lineal multivariante (modelo saturado)",
       subtitle = "Los valores predichos deberían estar cercanos a la diagonal",
       x = "Predicciones",
       y = "Valores reales")

g2 <- pred_test |> 
  select(.pred, SalePrice) |> 
  gather(Distribución, value) |> 
  ggplot(aes(x = value, color = Distribución, fill = Distribución)) + 
  geom_density(alpha = 0.6) + 
  theme_minimal() + 
  labs(title = "Distribución de las predicciones sobre los valores reales de SalePrice",
       x = "Distribuciones",
       y = "Frecuencia")

multiplot(g1, g2)

Como se puede apreciar en los gráficos, las predicciones son bastante buenas, bastante mejores que en el modelo univariante. A pesar de que se han incumplido algunos de los supuestos de regresión, la recta se ajusta bastante bien a los datos. A continuación, crearemos un modelo a caballo entre el univariante y el saturado, incluyendo únicamente en el modelo una selección del total de variables del dataset.

Fase 4.5: Receta, Modelo y Flujo (Regresión multivariante con selección de modelos)

Para este último modelo de regresión vamos a emplear únicamente un conjunto de variables seleccionadas a través de distintos criterios de información.

Seleccion de modelos

Para la selección de modelos, pasaremos directamente la receta a lm() para proceder a la regresión contra las variables escogidas.

house_prep <- 
  bake(house_rec_multi |> prep(), new_data = NULL)
ajuste_house_ulti <- 
  lm(data = house_prep , SalePrice ~ .)

En este caso, y tras varias pruebas, hemos decidido emplear el criterio de información BIC, por lo que le aplicaremos a la función stepAIC su penalización correspondiente:

set.seed(05492)

modBIC <-
  MASS::stepAIC(ajuste_house_ulti , k = log(nrow(house_train)))
Start:  AIC=-5503.97
SalePrice ~ LotArea + OverallCond + ExterCond + BsmtQual + BsmtCond + 
    HeatingQC + `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + 
    GarageCars + GarageQual + RemodAdd + New + Age + MSSubClass_X.150k + 
    MSZoning_RL + Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + 
    Neighborhood_X150k.200k + Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + 
    RoofStyle_Gam_Gable + Exterior1st_X.150k + CentralAir_Y + 
    Electrical_SBrkr + Functional_Typ + PavedDrive_Y + Fence_MnPrv + 
    MoSold_Feb_Aug_Sep_Nov_Dec + MoSold_Mar_Jun_Jul + YrSold_X2010 + 
    SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- MoSold_Mar_Jun_Jul          1   0.00004 28.267 -5511.2
- Exterior1st_X.150k          1   0.00011 28.267 -5511.2
- Electrical_SBrkr            1   0.00021 28.267 -5511.2
- ExterCond                   1   0.00142 28.268 -5511.2
- YrSold_X2010                1   0.00493 28.272 -5511.0
- New                         1   0.00625 28.273 -5510.9
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01073 28.277 -5510.7
- Fence_MnPrv                 1   0.01344 28.280 -5510.6
- BsmtCond                    1   0.02185 28.288 -5510.1
- LandSlope_Mod_Sev           1   0.03171 28.298 -5509.6
- MSSubClass_X.150k           1   0.03942 28.306 -5509.2
- `2ndFlrSF`                  1   0.05127 28.318 -5508.6
- Alley_None                  1   0.06139 28.328 -5508.1
- LandContour_Lvl             1   0.06187 28.328 -5508.1
- RoofStyle_Gam_Gable         1   0.08215 28.349 -5507.0
- `1stFlrSF`                  1   0.10321 28.370 -5505.9
<none>                                    28.267 -5504.0
- PavedDrive_Y                1   0.17295 28.440 -5502.3
- Neighborhood_X150k.200k     1   0.22259 28.489 -5499.8
- BldgType_fmCon_Dup_Twnhs    1   0.22948 28.496 -5499.4
- LotArea                     1   0.25478 28.521 -5498.2
- GarageQual                  1   0.27613 28.543 -5497.1
- CentralAir_Y                1   0.28574 28.552 -5496.6
- RemodAdd                    1   0.33524 28.602 -5494.0
- SaleCondition_Normal        1   0.33567 28.602 -5494.0
- HeatingQC                   1   0.40587 28.672 -5490.4
- MSZoning_RL                 1   0.52610 28.793 -5484.3
- FireplaceQu                 1   0.63386 28.901 -5478.9
- Age                         1   0.93476 29.201 -5463.8
- Functional_Typ              1   1.04377 29.310 -5458.3
- GarageCars                  1   1.57238 29.839 -5432.2
- GrLivArea                   1   1.93106 30.198 -5414.8
- OverallCond                 1   1.95981 30.227 -5413.4
- BsmtQual                    1   2.96959 31.236 -5365.4
- Neighborhood_X.200k         1   3.12340 31.390 -5358.2

Step:  AIC=-5511.25
SalePrice ~ LotArea + OverallCond + ExterCond + BsmtQual + BsmtCond + 
    HeatingQC + `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + 
    GarageCars + GarageQual + RemodAdd + New + Age + MSSubClass_X.150k + 
    MSZoning_RL + Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + 
    Neighborhood_X150k.200k + Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + 
    RoofStyle_Gam_Gable + Exterior1st_X.150k + CentralAir_Y + 
    Electrical_SBrkr + Functional_Typ + PavedDrive_Y + Fence_MnPrv + 
    MoSold_Feb_Aug_Sep_Nov_Dec + YrSold_X2010 + SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- Exterior1st_X.150k          1   0.00011 28.267 -5518.5
- Electrical_SBrkr            1   0.00021 28.267 -5518.5
- ExterCond                   1   0.00141 28.268 -5518.5
- YrSold_X2010                1   0.00511 28.272 -5518.3
- New                         1   0.00627 28.273 -5518.2
- Fence_MnPrv                 1   0.01342 28.280 -5517.8
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01343 28.280 -5517.8
- BsmtCond                    1   0.02188 28.289 -5517.4
- LandSlope_Mod_Sev           1   0.03197 28.299 -5516.9
- MSSubClass_X.150k           1   0.03938 28.306 -5516.5
- `2ndFlrSF`                  1   0.05123 28.318 -5515.9
- Alley_None                  1   0.06164 28.328 -5515.4
- LandContour_Lvl             1   0.06209 28.329 -5515.3
- RoofStyle_Gam_Gable         1   0.08215 28.349 -5514.3
- `1stFlrSF`                  1   0.10353 28.370 -5513.2
<none>                                    28.267 -5511.2
- PavedDrive_Y                1   0.17292 28.440 -5509.6
- Neighborhood_X150k.200k     1   0.22263 28.489 -5507.1
- BldgType_fmCon_Dup_Twnhs    1   0.22944 28.496 -5506.7
- LotArea                     1   0.25492 28.522 -5505.4
- GarageQual                  1   0.27609 28.543 -5504.3
- CentralAir_Y                1   0.28666 28.553 -5503.8
- RemodAdd                    1   0.33528 28.602 -5501.3
- SaleCondition_Normal        1   0.33564 28.602 -5501.3
- HeatingQC                   1   0.40645 28.673 -5497.7
- MSZoning_RL                 1   0.52662 28.793 -5491.6
- FireplaceQu                 1   0.63444 28.901 -5486.1
- Age                         1   0.93473 29.201 -5471.0
- Functional_Typ              1   1.04385 29.311 -5465.6
- GarageCars                  1   1.57360 29.840 -5439.4
- GrLivArea                   1   1.93324 30.200 -5421.9
- OverallCond                 1   1.96045 30.227 -5420.6
- BsmtQual                    1   2.98093 31.248 -5372.2
- Neighborhood_X.200k         1   3.15004 31.417 -5364.3

Step:  AIC=-5518.53
SalePrice ~ LotArea + OverallCond + ExterCond + BsmtQual + BsmtCond + 
    HeatingQC + `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + 
    GarageCars + GarageQual + RemodAdd + New + Age + MSSubClass_X.150k + 
    MSZoning_RL + Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + 
    Neighborhood_X150k.200k + Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + 
    RoofStyle_Gam_Gable + CentralAir_Y + Electrical_SBrkr + Functional_Typ + 
    PavedDrive_Y + Fence_MnPrv + MoSold_Feb_Aug_Sep_Nov_Dec + 
    YrSold_X2010 + SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- Electrical_SBrkr            1   0.00020 28.267 -5525.8
- ExterCond                   1   0.00142 28.268 -5525.7
- YrSold_X2010                1   0.00507 28.272 -5525.6
- New                         1   0.00640 28.273 -5525.5
- Fence_MnPrv                 1   0.01335 28.280 -5525.1
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01339 28.280 -5525.1
- BsmtCond                    1   0.02193 28.289 -5524.7
- LandSlope_Mod_Sev           1   0.03222 28.299 -5524.2
- MSSubClass_X.150k           1   0.03929 28.306 -5523.8
- `2ndFlrSF`                  1   0.05128 28.318 -5523.2
- Alley_None                  1   0.06154 28.328 -5522.6
- LandContour_Lvl             1   0.06199 28.329 -5522.6
- RoofStyle_Gam_Gable         1   0.08334 28.350 -5521.5
- `1stFlrSF`                  1   0.10349 28.370 -5520.5
<none>                                    28.267 -5518.5
- PavedDrive_Y                1   0.17313 28.440 -5516.9
- Neighborhood_X150k.200k     1   0.22541 28.492 -5514.2
- BldgType_fmCon_Dup_Twnhs    1   0.22971 28.497 -5514.0
- LotArea                     1   0.25492 28.522 -5512.7
- GarageQual                  1   0.27691 28.544 -5511.6
- CentralAir_Y                1   0.28673 28.553 -5511.1
- SaleCondition_Normal        1   0.33788 28.605 -5508.5
- RemodAdd                    1   0.33860 28.605 -5508.4
- HeatingQC                   1   0.40985 28.677 -5504.8
- MSZoning_RL                 1   0.52665 28.794 -5498.9
- FireplaceQu                 1   0.63457 28.901 -5493.4
- Age                         1   0.97330 29.240 -5476.4
- Functional_Typ              1   1.04383 29.311 -5472.9
- GarageCars                  1   1.58183 29.849 -5446.3
- GrLivArea                   1   1.93397 30.201 -5429.2
- OverallCond                 1   1.98051 30.247 -5426.9
- BsmtQual                    1   2.99256 31.259 -5378.9
- Neighborhood_X.200k         1   3.15151 31.418 -5371.5

Step:  AIC=-5525.81
SalePrice ~ LotArea + OverallCond + ExterCond + BsmtQual + BsmtCond + 
    HeatingQC + `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + 
    GarageCars + GarageQual + RemodAdd + New + Age + MSSubClass_X.150k + 
    MSZoning_RL + Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + 
    Neighborhood_X150k.200k + Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + 
    RoofStyle_Gam_Gable + CentralAir_Y + Functional_Typ + PavedDrive_Y + 
    Fence_MnPrv + MoSold_Feb_Aug_Sep_Nov_Dec + YrSold_X2010 + 
    SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- ExterCond                   1   0.00134 28.268 -5533.0
- YrSold_X2010                1   0.00516 28.272 -5532.8
- New                         1   0.00633 28.273 -5532.8
- Fence_MnPrv                 1   0.01320 28.280 -5532.4
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01345 28.280 -5532.4
- BsmtCond                    1   0.02218 28.289 -5531.9
- LandSlope_Mod_Sev           1   0.03235 28.299 -5531.4
- MSSubClass_X.150k           1   0.03930 28.306 -5531.1
- `2ndFlrSF`                  1   0.05148 28.319 -5530.4
- Alley_None                  1   0.06159 28.329 -5529.9
- LandContour_Lvl             1   0.06220 28.329 -5529.9
- RoofStyle_Gam_Gable         1   0.08345 28.351 -5528.8
- `1stFlrSF`                  1   0.10329 28.370 -5527.8
<none>                                    28.267 -5525.8
- PavedDrive_Y                1   0.17396 28.441 -5524.1
- Neighborhood_X150k.200k     1   0.22748 28.494 -5521.4
- BldgType_fmCon_Dup_Twnhs    1   0.22963 28.497 -5521.3
- LotArea                     1   0.25472 28.522 -5520.0
- GarageQual                  1   0.27719 28.544 -5518.8
- CentralAir_Y                1   0.29888 28.566 -5517.7
- SaleCondition_Normal        1   0.33835 28.605 -5515.7
- RemodAdd                    1   0.34266 28.610 -5515.5
- HeatingQC                   1   0.40982 28.677 -5512.1
- MSZoning_RL                 1   0.52933 28.796 -5506.0
- FireplaceQu                 1   0.63685 28.904 -5500.6
- Age                         1   0.99837 29.265 -5482.4
- Functional_Typ              1   1.04408 29.311 -5480.1
- GarageCars                  1   1.58615 29.853 -5453.4
- GrLivArea                   1   1.93726 30.204 -5436.3
- OverallCond                 1   1.99445 30.262 -5433.6
- BsmtQual                    1   3.00017 31.267 -5385.8
- Neighborhood_X.200k         1   3.15292 31.420 -5378.7

Step:  AIC=-5533.02
SalePrice ~ LotArea + OverallCond + BsmtQual + BsmtCond + HeatingQC + 
    `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + 
    GarageQual + RemodAdd + New + Age + MSSubClass_X.150k + MSZoning_RL + 
    Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + Fence_MnPrv + 
    MoSold_Feb_Aug_Sep_Nov_Dec + YrSold_X2010 + SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- YrSold_X2010                1   0.00516 28.273 -5540.0
- New                         1   0.00663 28.275 -5540.0
- Fence_MnPrv                 1   0.01288 28.281 -5539.6
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01359 28.282 -5539.6
- BsmtCond                    1   0.02115 28.290 -5539.2
- LandSlope_Mod_Sev           1   0.03294 28.301 -5538.6
- MSSubClass_X.150k           1   0.03906 28.307 -5538.3
- `2ndFlrSF`                  1   0.05255 28.321 -5537.6
- Alley_None                  1   0.06114 28.329 -5537.2
- LandContour_Lvl             1   0.06223 28.331 -5537.1
- RoofStyle_Gam_Gable         1   0.08356 28.352 -5536.0
- `1stFlrSF`                  1   0.10220 28.371 -5535.0
<none>                                    28.268 -5533.0
- PavedDrive_Y                1   0.17270 28.441 -5531.4
- Neighborhood_X150k.200k     1   0.22782 28.496 -5528.6
- BldgType_fmCon_Dup_Twnhs    1   0.22868 28.497 -5528.5
- LotArea                     1   0.25391 28.522 -5527.3
- GarageQual                  1   0.27594 28.544 -5526.1
- CentralAir_Y                1   0.29946 28.568 -5524.9
- SaleCondition_Normal        1   0.33894 28.607 -5522.9
- RemodAdd                    1   0.34162 28.610 -5522.8
- HeatingQC                   1   0.40887 28.677 -5519.3
- MSZoning_RL                 1   0.53406 28.802 -5513.0
- FireplaceQu                 1   0.63748 28.906 -5507.8
- Age                         1   0.99715 29.265 -5489.7
- Functional_Typ              1   1.04323 29.312 -5487.4
- GarageCars                  1   1.58557 29.854 -5460.6
- GrLivArea                   1   1.94693 30.215 -5443.1
- OverallCond                 1   2.17067 30.439 -5432.3
- BsmtQual                    1   3.01602 31.284 -5392.3
- Neighborhood_X.200k         1   3.15315 31.422 -5385.9

Step:  AIC=-5540.04
SalePrice ~ LotArea + OverallCond + BsmtQual + BsmtCond + HeatingQC + 
    `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + 
    GarageQual + RemodAdd + New + Age + MSSubClass_X.150k + MSZoning_RL + 
    Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + Fence_MnPrv + 
    MoSold_Feb_Aug_Sep_Nov_Dec + SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- New                         1   0.00637 28.280 -5547.0
- Fence_MnPrv                 1   0.01263 28.286 -5546.7
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01625 28.290 -5546.5
- BsmtCond                    1   0.02068 28.294 -5546.3
- LandSlope_Mod_Sev           1   0.03332 28.307 -5545.6
- MSSubClass_X.150k           1   0.03920 28.313 -5545.3
- `2ndFlrSF`                  1   0.05331 28.327 -5544.6
- Alley_None                  1   0.06237 28.336 -5544.1
- LandContour_Lvl             1   0.06342 28.337 -5544.1
- RoofStyle_Gam_Gable         1   0.08270 28.356 -5543.1
- `1stFlrSF`                  1   0.10160 28.375 -5542.1
<none>                                    28.273 -5540.0
- PavedDrive_Y                1   0.17423 28.448 -5538.4
- Neighborhood_X150k.200k     1   0.22636 28.500 -5535.7
- BldgType_fmCon_Dup_Twnhs    1   0.22813 28.502 -5535.6
- LotArea                     1   0.25428 28.528 -5534.3
- GarageQual                  1   0.27578 28.549 -5533.2
- CentralAir_Y                1   0.30002 28.573 -5531.9
- SaleCondition_Normal        1   0.34136 28.615 -5529.8
- RemodAdd                    1   0.34178 28.615 -5529.8
- HeatingQC                   1   0.40673 28.680 -5526.5
- MSZoning_RL                 1   0.53554 28.809 -5519.9
- FireplaceQu                 1   0.63340 28.907 -5515.0
- Age                         1   0.99472 29.268 -5496.8
- Functional_Typ              1   1.04195 29.316 -5494.5
- GarageCars                  1   1.58222 29.856 -5467.8
- GrLivArea                   1   1.95301 30.227 -5449.8
- OverallCond                 1   2.17522 30.449 -5439.1
- BsmtQual                    1   3.02848 31.302 -5398.8
- Neighborhood_X.200k         1   3.15576 31.429 -5392.8

Step:  AIC=-5547
SalePrice ~ LotArea + OverallCond + BsmtQual + BsmtCond + HeatingQC + 
    `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + 
    GarageQual + RemodAdd + Age + MSSubClass_X.150k + MSZoning_RL + 
    Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + Fence_MnPrv + 
    MoSold_Feb_Aug_Sep_Nov_Dec + SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- Fence_MnPrv                 1   0.01289 28.293 -5553.6
- MoSold_Feb_Aug_Sep_Nov_Dec  1   0.01468 28.295 -5553.5
- BsmtCond                    1   0.02236 28.302 -5553.1
- LandSlope_Mod_Sev           1   0.03181 28.312 -5552.6
- MSSubClass_X.150k           1   0.03872 28.319 -5552.3
- `2ndFlrSF`                  1   0.05226 28.332 -5551.6
- LandContour_Lvl             1   0.06144 28.341 -5551.1
- Alley_None                  1   0.06193 28.342 -5551.1
- RoofStyle_Gam_Gable         1   0.08418 28.364 -5549.9
- `1stFlrSF`                  1   0.10241 28.382 -5549.0
<none>                                    28.280 -5547.0
- PavedDrive_Y                1   0.17320 28.453 -5545.4
- Neighborhood_X150k.200k     1   0.22173 28.502 -5542.9
- BldgType_fmCon_Dup_Twnhs    1   0.23447 28.514 -5542.2
- LotArea                     1   0.26073 28.541 -5540.9
- GarageQual                  1   0.27658 28.556 -5540.1
- CentralAir_Y                1   0.29528 28.575 -5539.1
- SaleCondition_Normal        1   0.35858 28.639 -5535.9
- RemodAdd                    1   0.40182 28.682 -5533.7
- HeatingQC                   1   0.40539 28.685 -5533.5
- MSZoning_RL                 1   0.53381 28.814 -5527.0
- FireplaceQu                 1   0.63335 28.913 -5521.9
- Functional_Typ              1   1.03985 29.320 -5501.6
- Age                         1   1.12845 29.408 -5497.2
- GarageCars                  1   1.57672 29.857 -5475.1
- GrLivArea                   1   1.94981 30.230 -5456.9
- OverallCond                 1   2.17327 30.453 -5446.2
- BsmtQual                    1   3.02553 31.305 -5405.9
- Neighborhood_X.200k         1   3.15298 31.433 -5400.0

Step:  AIC=-5553.62
SalePrice ~ LotArea + OverallCond + BsmtQual + BsmtCond + HeatingQC + 
    `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + 
    GarageQual + RemodAdd + Age + MSSubClass_X.150k + MSZoning_RL + 
    Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + MoSold_Feb_Aug_Sep_Nov_Dec + 
    SaleCondition_Normal

                             Df Sum of Sq    RSS     AIC
- MoSold_Feb_Aug_Sep_Nov_Dec  1    0.0155 28.308 -5560.1
- BsmtCond                    1    0.0211 28.314 -5559.8
- LandSlope_Mod_Sev           1    0.0321 28.325 -5559.2
- MSSubClass_X.150k           1    0.0391 28.332 -5558.9
- `2ndFlrSF`                  1    0.0529 28.346 -5558.2
- LandContour_Lvl             1    0.0588 28.352 -5557.9
- Alley_None                  1    0.0618 28.355 -5557.7
- RoofStyle_Gam_Gable         1    0.0829 28.376 -5556.6
- `1stFlrSF`                  1    0.1007 28.393 -5555.7
<none>                                    28.293 -5553.6
- PavedDrive_Y                1    0.1715 28.464 -5552.1
- Neighborhood_X150k.200k     1    0.2220 28.515 -5549.5
- BldgType_fmCon_Dup_Twnhs    1    0.2257 28.518 -5549.3
- LotArea                     1    0.2610 28.554 -5547.5
- GarageQual                  1    0.2770 28.570 -5546.7
- CentralAir_Y                1    0.2904 28.583 -5546.0
- SaleCondition_Normal        1    0.3647 28.657 -5542.2
- RemodAdd                    1    0.4094 28.702 -5539.9
- HeatingQC                   1    0.4264 28.719 -5539.1
- MSZoning_RL                 1    0.5292 28.822 -5533.8
- FireplaceQu                 1    0.6418 28.934 -5528.2
- Functional_Typ              1    1.0614 29.354 -5507.1
- Age                         1    1.1544 29.447 -5502.5
- GarageCars                  1    1.5925 29.885 -5481.0
- GrLivArea                   1    1.9596 30.252 -5463.1
- OverallCond                 1    2.1604 30.453 -5453.5
- BsmtQual                    1    3.0328 31.326 -5412.2
- Neighborhood_X.200k         1    3.1684 31.461 -5405.9

Step:  AIC=-5560.11
SalePrice ~ LotArea + OverallCond + BsmtQual + BsmtCond + HeatingQC + 
    `1stFlrSF` + `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + 
    GarageQual + RemodAdd + Age + MSSubClass_X.150k + MSZoning_RL + 
    Alley_None + LandContour_Lvl + LandSlope_Mod_Sev + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- BsmtCond                  1   0.01902 28.327 -5566.4
- LandSlope_Mod_Sev         1   0.03193 28.340 -5565.7
- MSSubClass_X.150k         1   0.03729 28.346 -5565.5
- `2ndFlrSF`                1   0.05455 28.363 -5564.6
- LandContour_Lvl           1   0.06185 28.370 -5564.2
- Alley_None                1   0.06240 28.371 -5564.2
- RoofStyle_Gam_Gable       1   0.08440 28.393 -5563.0
- `1stFlrSF`                1   0.09833 28.407 -5562.3
<none>                                  28.308 -5560.1
- PavedDrive_Y              1   0.17043 28.479 -5558.6
- Neighborhood_X150k.200k   1   0.21608 28.524 -5556.3
- BldgType_fmCon_Dup_Twnhs  1   0.22927 28.538 -5555.6
- LotArea                   1   0.26051 28.569 -5554.0
- GarageQual                1   0.27898 28.587 -5553.1
- CentralAir_Y              1   0.29039 28.599 -5552.5
- SaleCondition_Normal      1   0.37453 28.683 -5548.2
- RemodAdd                  1   0.40543 28.714 -5546.6
- HeatingQC                 1   0.43295 28.741 -5545.2
- MSZoning_RL               1   0.53667 28.845 -5540.0
- FireplaceQu               1   0.63838 28.947 -5534.8
- Functional_Typ            1   1.06669 29.375 -5513.4
- Age                       1   1.14907 29.457 -5509.3
- GarageCars                1   1.60586 29.914 -5486.8
- GrLivArea                 1   1.96752 30.276 -5469.3
- OverallCond               1   2.17026 30.479 -5459.5
- BsmtQual                  1   3.03044 31.339 -5418.9
- Neighborhood_X.200k       1   3.15411 31.462 -5413.2

Step:  AIC=-5566.41
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `1stFlrSF` + 
    `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + GarageQual + 
    RemodAdd + Age + MSSubClass_X.150k + MSZoning_RL + Alley_None + 
    LandContour_Lvl + LandSlope_Mod_Sev + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- MSSubClass_X.150k         1    0.0358 28.363 -5571.9
- LandSlope_Mod_Sev         1    0.0363 28.364 -5571.8
- `2ndFlrSF`                1    0.0554 28.383 -5570.8
- Alley_None                1    0.0598 28.387 -5570.6
- LandContour_Lvl           1    0.0661 28.393 -5570.3
- RoofStyle_Gam_Gable       1    0.0829 28.410 -5569.4
- `1stFlrSF`                1    0.0963 28.424 -5568.7
<none>                                  28.327 -5566.4
- PavedDrive_Y              1    0.1769 28.504 -5564.6
- Neighborhood_X150k.200k   1    0.2105 28.538 -5562.9
- BldgType_fmCon_Dup_Twnhs  1    0.2344 28.562 -5561.7
- LotArea                   1    0.2648 28.592 -5560.1
- GarageQual                1    0.2789 28.606 -5559.4
- CentralAir_Y              1    0.3085 28.636 -5557.9
- SaleCondition_Normal      1    0.3654 28.693 -5555.0
- RemodAdd                  1    0.4118 28.739 -5552.6
- HeatingQC                 1    0.4320 28.759 -5551.6
- MSZoning_RL               1    0.5416 28.869 -5546.1
- FireplaceQu               1    0.6371 28.964 -5541.2
- Functional_Typ            1    1.0810 29.408 -5519.0
- Age                       1    1.1342 29.461 -5516.4
- GarageCars                1    1.6032 29.930 -5493.3
- GrLivArea                 1    1.9723 30.300 -5475.4
- OverallCond               1    2.3591 30.686 -5456.9
- Neighborhood_X.200k       1    3.1352 31.462 -5420.4
- BsmtQual                  1    3.6434 31.971 -5397.0

Step:  AIC=-5571.86
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `1stFlrSF` + 
    `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + GarageQual + 
    RemodAdd + Age + MSZoning_RL + Alley_None + LandContour_Lvl + 
    LandSlope_Mod_Sev + Neighborhood_X150k.200k + Neighborhood_X.200k + 
    BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + CentralAir_Y + 
    Functional_Typ + PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- LandSlope_Mod_Sev         1    0.0348 28.398 -5577.4
- Alley_None                1    0.0675 28.431 -5575.7
- LandContour_Lvl           1    0.0690 28.432 -5575.6
- `1stFlrSF`                1    0.0761 28.439 -5575.2
- `2ndFlrSF`                1    0.0830 28.446 -5574.9
- RoofStyle_Gam_Gable       1    0.0917 28.455 -5574.4
<none>                                  28.363 -5571.9
- PavedDrive_Y              1    0.1834 28.547 -5569.7
- Neighborhood_X150k.200k   1    0.2250 28.588 -5567.6
- BldgType_fmCon_Dup_Twnhs  1    0.2544 28.617 -5566.1
- GarageQual                1    0.2768 28.640 -5565.0
- CentralAir_Y              1    0.3020 28.665 -5563.7
- LotArea                   1    0.3418 28.705 -5561.7
- SaleCondition_Normal      1    0.3601 28.723 -5560.7
- HeatingQC                 1    0.4229 28.786 -5557.5
- RemodAdd                  1    0.4519 28.815 -5556.1
- MSZoning_RL               1    0.5966 28.960 -5548.8
- FireplaceQu               1    0.6122 28.975 -5548.0
- Functional_Typ            1    1.1342 29.497 -5521.9
- Age                       1    1.1918 29.555 -5519.0
- GarageCars                1    1.6328 29.996 -5497.4
- GrLivArea                 1    2.2809 30.644 -5466.2
- OverallCond               1    2.3719 30.735 -5461.9
- Neighborhood_X.200k       1    3.1067 31.470 -5427.4
- BsmtQual                  1    3.6077 31.971 -5404.3

Step:  AIC=-5577.35
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `1stFlrSF` + 
    `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + GarageQual + 
    RemodAdd + Age + MSZoning_RL + Alley_None + LandContour_Lvl + 
    Neighborhood_X150k.200k + Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + 
    RoofStyle_Gam_Gable + CentralAir_Y + Functional_Typ + PavedDrive_Y + 
    SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- LandContour_Lvl           1    0.0383 28.436 -5582.7
- Alley_None                1    0.0709 28.469 -5581.0
- `1stFlrSF`                1    0.0757 28.474 -5580.8
- `2ndFlrSF`                1    0.0828 28.481 -5580.4
- RoofStyle_Gam_Gable       1    0.0917 28.490 -5579.9
<none>                                  28.398 -5577.4
- PavedDrive_Y              1    0.1894 28.587 -5574.9
- Neighborhood_X150k.200k   1    0.2311 28.629 -5572.8
- BldgType_fmCon_Dup_Twnhs  1    0.2540 28.652 -5571.6
- GarageQual                1    0.2787 28.677 -5570.4
- CentralAir_Y              1    0.3083 28.706 -5568.9
- SaleCondition_Normal      1    0.3638 28.762 -5566.1
- LotArea                   1    0.3738 28.772 -5565.5
- HeatingQC                 1    0.4211 28.819 -5563.1
- RemodAdd                  1    0.4524 28.850 -5561.6
- MSZoning_RL               1    0.5869 28.985 -5554.8
- FireplaceQu               1    0.6097 29.008 -5553.6
- Functional_Typ            1    1.1050 29.503 -5528.9
- Age                       1    1.1746 29.572 -5525.5
- GarageCars                1    1.6234 30.021 -5503.5
- GrLivArea                 1    2.2737 30.672 -5472.2
- OverallCond               1    2.3779 30.776 -5467.2
- Neighborhood_X.200k       1    3.1839 31.582 -5429.5
- BsmtQual                  1    3.6144 32.012 -5409.7

Step:  AIC=-5582.67
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `1stFlrSF` + 
    `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + GarageQual + 
    RemodAdd + Age + MSZoning_RL + Alley_None + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- `1stFlrSF`                1    0.0726 28.509 -5586.2
- Alley_None                1    0.0751 28.511 -5586.1
- `2ndFlrSF`                1    0.0868 28.523 -5585.5
- RoofStyle_Gam_Gable       1    0.0906 28.527 -5585.3
<none>                                  28.436 -5582.7
- PavedDrive_Y              1    0.2008 28.637 -5579.7
- Neighborhood_X150k.200k   1    0.2206 28.657 -5578.7
- BldgType_fmCon_Dup_Twnhs  1    0.2683 28.705 -5576.2
- GarageQual                1    0.2833 28.720 -5575.5
- CentralAir_Y              1    0.3107 28.747 -5574.1
- LotArea                   1    0.3426 28.779 -5572.5
- SaleCondition_Normal      1    0.3732 28.809 -5570.9
- HeatingQC                 1    0.4342 28.870 -5567.8
- RemodAdd                  1    0.4627 28.899 -5566.4
- FireplaceQu               1    0.5972 29.033 -5559.6
- MSZoning_RL               1    0.6179 29.054 -5558.6
- Functional_Typ            1    1.1186 29.555 -5533.6
- Age                       1    1.1928 29.629 -5530.0
- GarageCars                1    1.6344 30.071 -5508.4
- GrLivArea                 1    2.3120 30.748 -5475.8
- OverallCond               1    2.3643 30.801 -5473.3
- Neighborhood_X.200k       1    3.1516 31.588 -5436.5
- BsmtQual                  1    3.6052 32.041 -5415.7

Step:  AIC=-5586.23
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `2ndFlrSF` + 
    GrLivArea + FireplaceQu + GarageCars + GarageQual + RemodAdd + 
    Age + MSZoning_RL + Alley_None + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + 
    CentralAir_Y + Functional_Typ + PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- Alley_None                1    0.0764 28.585 -5589.6
- RoofStyle_Gam_Gable       1    0.0887 28.598 -5589.0
<none>                                  28.509 -5586.2
- PavedDrive_Y              1    0.1960 28.705 -5583.5
- Neighborhood_X150k.200k   1    0.1999 28.709 -5583.3
- BldgType_fmCon_Dup_Twnhs  1    0.2748 28.784 -5579.5
- GarageQual                1    0.2902 28.799 -5578.7
- CentralAir_Y              1    0.3214 28.830 -5577.1
- SaleCondition_Normal      1    0.3624 28.871 -5575.1
- LotArea                   1    0.3915 28.900 -5573.6
- HeatingQC                 1    0.4295 28.938 -5571.7
- RemodAdd                  1    0.4403 28.949 -5571.1
- FireplaceQu               1    0.6183 29.127 -5562.2
- MSZoning_RL               1    0.6194 29.128 -5562.1
- Functional_Typ            1    1.0849 29.594 -5539.0
- Age                       1    1.1662 29.675 -5535.0
- GarageCars                1    1.6322 30.141 -5512.2
- `2ndFlrSF`                1    2.0406 30.549 -5492.6
- OverallCond               1    2.4071 30.916 -5475.2
- Neighborhood_X.200k       1    3.1062 31.615 -5442.5
- BsmtQual                  1    3.6271 32.136 -5418.7
- GrLivArea                 1   17.3706 45.879 -4898.8

Step:  AIC=-5589.61
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `2ndFlrSF` + 
    GrLivArea + FireplaceQu + GarageCars + GarageQual + RemodAdd + 
    Age + MSZoning_RL + Neighborhood_X150k.200k + Neighborhood_X.200k + 
    BldgType_fmCon_Dup_Twnhs + RoofStyle_Gam_Gable + CentralAir_Y + 
    Functional_Typ + PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
- RoofStyle_Gam_Gable       1    0.0959 28.681 -5592.0
<none>                                  28.585 -5589.6
- Neighborhood_X150k.200k   1    0.2151 28.800 -5586.0
- PavedDrive_Y              1    0.2319 28.817 -5585.1
- BldgType_fmCon_Dup_Twnhs  1    0.2686 28.854 -5583.2
- GarageQual                1    0.2872 28.872 -5582.3
- CentralAir_Y              1    0.3544 28.940 -5578.9
- SaleCondition_Normal      1    0.3742 28.959 -5577.9
- HeatingQC                 1    0.4151 29.000 -5575.8
- LotArea                   1    0.4220 29.007 -5575.5
- RemodAdd                  1    0.4618 29.047 -5573.5
- FireplaceQu               1    0.6560 29.241 -5563.8
- MSZoning_RL               1    0.6669 29.252 -5563.2
- Functional_Typ            1    1.0892 29.674 -5542.3
- Age                       1    1.1631 29.748 -5538.7
- GarageCars                1    1.6045 30.190 -5517.2
- `2ndFlrSF`                1    2.1151 30.700 -5492.7
- OverallCond               1    2.3804 30.966 -5480.1
- Neighborhood_X.200k       1    3.0737 31.659 -5447.8
- BsmtQual                  1    3.6742 32.259 -5420.4
- GrLivArea                 1   17.3370 45.922 -4904.8

Step:  AIC=-5592.01
SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + `2ndFlrSF` + 
    GrLivArea + FireplaceQu + GarageCars + GarageQual + RemodAdd + 
    Age + MSZoning_RL + Neighborhood_X150k.200k + Neighborhood_X.200k + 
    BldgType_fmCon_Dup_Twnhs + CentralAir_Y + Functional_Typ + 
    PavedDrive_Y + SaleCondition_Normal

                           Df Sum of Sq    RSS     AIC
<none>                                  28.681 -5592.0
- Neighborhood_X150k.200k   1    0.1852 28.866 -5589.9
- PavedDrive_Y              1    0.2269 28.908 -5587.8
- BldgType_fmCon_Dup_Twnhs  1    0.2758 28.957 -5585.3
- GarageQual                1    0.2927 28.974 -5584.5
- CentralAir_Y              1    0.3404 29.021 -5582.1
- SaleCondition_Normal      1    0.3796 29.061 -5580.1
- HeatingQC                 1    0.4104 29.092 -5578.6
- LotArea                   1    0.4502 29.131 -5576.6
- RemodAdd                  1    0.4972 29.178 -5574.2
- MSZoning_RL               1    0.6623 29.343 -5566.0
- FireplaceQu               1    0.6813 29.362 -5565.0
- Functional_Typ            1    1.0742 29.755 -5545.6
- Age                       1    1.1514 29.832 -5541.8
- GarageCars                1    1.6213 30.302 -5519.0
- OverallCond               1    2.4337 31.115 -5480.4
- `2ndFlrSF`                1    2.4685 31.150 -5478.8
- Neighborhood_X.200k       1    3.0523 31.733 -5451.6
- BsmtQual                  1    3.7266 32.408 -5420.9
- GrLivArea                 1   18.2085 46.890 -4881.6

Resumimos el modelo tras su procesado con summary():

summary(modBIC)

Call:
lm(formula = SalePrice ~ LotArea + OverallCond + BsmtQual + HeatingQC + 
    `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + GarageQual + 
    RemodAdd + Age + MSZoning_RL + Neighborhood_X150k.200k + 
    Neighborhood_X.200k + BldgType_fmCon_Dup_Twnhs + CentralAir_Y + 
    Functional_Typ + PavedDrive_Y + SaleCondition_Normal, data = house_prep)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.13018 -0.07507  0.00400  0.08492  0.50303 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)              11.663334   0.026693 436.937  < 2e-16 ***
LotArea                   0.022914   0.004820   4.754 2.19e-06 ***
OverallCond               0.047640   0.004310  11.054  < 2e-16 ***
BsmtQual                  0.074346   0.005435  13.679  < 2e-16 ***
HeatingQC                 0.021426   0.004720   4.539 6.11e-06 ***
`2ndFlrSF`               -0.057006   0.005121 -11.133  < 2e-16 ***
GrLivArea                 0.201236   0.006656  30.236  < 2e-16 ***
FireplaceQu               0.026632   0.004554   5.849 6.12e-09 ***
GarageCars                0.049482   0.005485   9.022  < 2e-16 ***
GarageQual                0.016549   0.004317   3.834 0.000132 ***
RemodAdd                 -0.021793   0.004362  -4.996 6.56e-07 ***
Age                      -0.042503   0.005590  -7.603 5.18e-14 ***
MSZoning_RL               0.070250   0.012182   5.766 9.89e-09 ***
Neighborhood_X150k.200k   0.033704   0.011054   3.049 0.002337 ** 
Neighborhood_X.200k       0.157659   0.012736  12.379  < 2e-16 ***
BldgType_fmCon_Dup_Twnhs -0.054103   0.014539  -3.721 0.000206 ***
CentralAir_Y              0.072612   0.017565   4.134 3.77e-05 ***
Functional_Typ            0.114182   0.015548   7.344 3.46e-13 ***
PavedDrive_Y              0.052063   0.015426   3.375 0.000758 ***
SaleCondition_Normal      0.044835   0.010270   4.365 1.36e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1411 on 1440 degrees of freedom
Multiple R-squared:  0.8768,    Adjusted R-squared:  0.8752 
F-statistic: 539.4 on 19 and 1440 DF,  p-value: < 2.2e-16

Como se puede observar, ha eliminado una gran cantidad de variables. Pasamos ahora a reescribir la receta únicamente con las variables seleccionadas:

Receta para el método de regresión multivariante con selección de modelos

Modificación previa de las variables según BIC

house_bic <- 
  house_train |> mutate(MSZoning = MSZoning == "RL")
house_bic <- 
  house_bic |> mutate(Neighborhood = Neighborhood == "X.200k" )
house_bic <- 
  house_bic |> mutate(BldgType = BldgType == "fmCon_Dup_Twnhs" )
house_bic <- 
  house_bic |> mutate(CentralAir = CentralAir == "Y")
house_bic <- 
  house_bic |> mutate(Functional = Functional == "Typ")
house_bic <- 
  house_bic |> mutate(PavedDrive = PavedDrive == "Y")
house_bic <- 
  house_bic |> mutate(SaleCondition = SaleCondition == "Normal")

Aplicación de roles

Volvemos a definir una nueva receta indicándole el conjunto donde tenemos validación y train, y enfrentando nuestra variable objetivo SalePrice a las variables seleccionadas por step_AIC(). Después, asignamos posibles roles, sujetos a modificación, que nos permitan diferenciar acciones entre las variables (sobre todo en la sección outliers).

# Receta
house_rec_multi_bic <-
  # Fórmula y datos
  recipe(data = house_bic, SalePrice ~ LotArea + OverallCond + BsmtQual + 
           HeatingQC + `2ndFlrSF` + GrLivArea + FireplaceQu + GarageCars + 
           GarageQual + RemodAdd + Age + MSZoning + Neighborhood + BldgType + 
           CentralAir + Functional + PavedDrive + SaleCondition)|>
  # Roles
  add_role(where(is.factor), 
           new_role = "cualitativa") |> 
  add_role(where(is.numeric), 
           new_role = "cuantitativa") |> 
  add_role(LotArea, `2ndFlrSF`, 
           new_role = "mediana") |> 
  add_role(all_nominal_predictors(), 
           new_role = "moda") |> 
  add_role(OverallCond,
           new_role = "imputar media")

Reagrupación de las variables cualitativas

En este modelo no será necesario reagrupar. El criterio de información BIC ya ha seleccionado las categorías más apropiadas, si las hubiera, de cada variable.

Recategorización de las variables cuantitativas

En los métodos de regresión, la recategorización de las variables cuantitativas no es estrictamente necesaria. Por tanto, en este receta, del total de variables cuantitativas no vamos a modificar ninguna.

Outliers

multiplot(box1, box2, box3, box4, box5, box6, box7, box8, box9, box10, cols = 2) 
multiplot(box11, box12, box13, box14, box15, box16, box17, cols = 2)

Si observamos estos gráficos de cajas y bigotes, todas nuestras variables cuantitativas continuas son asimétricas, por lo que se detectarán los outliers y se imputarán los ausentes por la mediana. Para el caso de las semi-cuantitativas (KitchenQual, por ejemplo) se imputarán directamente los ausentes por uno de los dos métodos en función de la simetría o asimetría de la variable. Para el resto de variables cualitativas, imputamos directamente por la moda.

house_rec_multi_bic <-
  house_rec_multi_bic |> 
  # Detección de outliers por la mediana y por la media
  step_mutate(across(has_role("mediana"), function(x) { ifelse(abs(scores(x, type = "mad")) > 3 & !is.na(x), NA, x) })) |>
  # Imputación de ausentes por la mediana, la media y la moda
  step_impute_median(has_role("mediana")) |>
  step_impute_mean(has_role("imputar_media")) |>
  step_impute_mode(has_role("moda")) |> 
  step_impute_knn(all_numeric_predictors() , -has_role("mediana"), -has_role("imputar_media") )

YeoJohnson para la normalidad de los residuos

A fin de que nuestro modelo cumpla los test de normalidad de los residuos, vamos a incluir en la receta la función step_YeoJohnson, que se emplea para eliminar la asimetría en los datos numéricos del modelo.

house_rec_multi_bic <-
  house_rec_multi_bic |> 
  step_YeoJohnson(all_numeric_predictors())

Transformación logarítmica de algunas variables

Transformamos algunas variables a logarítmicas para mejorar su relación.

house_rec_multi_bic <-
  house_rec_multi_bic |> 
  step_log(SalePrice, offset = 1)

Normalizar por rango

Normalizamos nuestras variables por rango para que todas tengan el mismo peso, entre 0 y 1.

house_rec_multi_bic <-
  house_rec_multi_bic |> 
  step_normalize(all_numeric_predictors())

Variables dummy

Como esta receta es para un modelo de regresión multivariante, debemos dummyficar nuestras variables cualitativas. Para ello, tomamos todas las nominales, menos nuestra variable objetivo.

house_rec_multi_bic <-
  house_rec_multi_bic |>
  step_dummy(all_nominal(), -all_outcomes())

Filtro de cero varianza

house_rec_multi_bic <-
  house_rec_multi_bic |>
  step_zv(all_predictors())

Filtro de correlación

Aplicamos el filtro de correlación a nuestras variables numéricas.

house_rec_multi_bic <-
  house_rec_multi_bic |> 
  step_corr(all_numeric_predictors(), threshold = 0.6)

Horneado

Por último, horneamos nuestra receta para comprobar que todas nuestras nuevas variables recategorizadas se hayan creado correctamente.

bake(house_rec_multi_bic |>  prep(), new_data = NULL)
# A tibble: 1,460 × 17
   LotArea OverallCond BsmtQual HeatingQC `2ndFlrSF` GrLivArea
     <dbl>       <dbl>    <dbl>     <dbl>      <dbl>     <dbl>
 1 -0.141       -0.477    0.576     0.930      1.17      0.528
 2  0.106        2.01     0.576     0.930     -0.871    -0.383
 3  0.414       -0.477    0.576     0.930      1.17      0.659
 4  0.0955      -0.477   -0.695    -0.330      1.15      0.541
 5  0.877       -0.477    0.576     0.930      1.21      1.28 
 6  0.857       -0.477    0.576     0.930      1.09     -0.154
 7  0.201       -0.477    2.13      0.930     -0.871     0.500
 8  0.257        0.440    0.576     0.930      1.20      1.13 
 9 -0.761       -0.477   -0.695    -0.330      1.15      0.639
10 -0.391        0.440   -0.695     0.930     -0.871    -0.857
# … with 1,450 more rows, and 11 more variables: FireplaceQu <dbl>,
#   GarageCars <dbl>, GarageQual <dbl>, RemodAdd <dbl>, Age <dbl>,
#   MSZoning <lgl>, CentralAir <lgl>, Functional <lgl>,
#   PavedDrive <lgl>, SaleCondition <lgl>, SalePrice <dbl>

Flujo y ajuste

# Creación del flujo
reg_wflow_bic <-
  workflow() |> 
  add_model(reg_lineal) |> 
  add_recipe(house_rec_multi_bic)
set.seed(05492)

reg_fit_bic <-
  reg_wflow_bic |> fit(data = house_train)
reg_fit_bic 
══ Workflow [trained] ════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ──────────────────────────────────────────────────────
11 Recipe Steps

• step_mutate()
• step_impute_median()
• step_impute_mean()
• step_impute_mode()
• step_impute_knn()
• step_YeoJohnson()
• step_log()
• step_normalize()
• step_dummy()
• step_zv()
• ...
• and 1 more step.

── Model ─────────────────────────────────────────────────────────────

Call:
stats::lm(formula = ..y ~ ., data = data)

Coefficients:
          (Intercept)                LotArea            OverallCond  
           11.8486103              0.0354377              0.0491536  
             BsmtQual              HeatingQC             `2ndFlrSF`  
            0.0673094              0.0183775             -0.0548430  
            GrLivArea            FireplaceQu             GarageCars  
            0.1960803              0.0267381              0.0432562  
           GarageQual               RemodAdd                    Age  
            0.0157258             -0.0167712             -0.0378894  
          MSZoning_FV            MSZoning_RH            MSZoning_RM  
            0.0797640              0.0090258              0.0296974  
 Neighborhood_Blueste    Neighborhood_BrDale   Neighborhood_BrkSide  
           -0.1233748             -0.0445626             -0.1180328  
 Neighborhood_ClearCr   Neighborhood_CollgCr   Neighborhood_Crawfor  
            0.0032435             -0.0258185              0.0083973  
 Neighborhood_Edwards   Neighborhood_Gilbert    Neighborhood_IDOTRR  
           -0.1625234             -0.0772346             -0.2696271  
 Neighborhood_MeadowV   Neighborhood_Mitchel     Neighborhood_NAmes  
           -0.1879102             -0.0919997             -0.0930020  
 Neighborhood_NoRidge   Neighborhood_NPkVill   Neighborhood_NridgHt  
            0.1261509             -0.0229875              0.1328166  
  Neighborhood_NWAmes   Neighborhood_OldTown    Neighborhood_Sawyer  
           -0.1066500             -0.2251197             -0.1199513  
 Neighborhood_SawyerW   Neighborhood_StoneBr     Neighborhood_SWISU  
           -0.0555471              0.2130249             -0.1420178  
  Neighborhood_Timber   Neighborhood_Veenker       BldgType_X2fmCon  
           -0.0080007              0.0485150             -0.0215432  
      BldgType_Duplex         BldgType_Twnhs        BldgType_TwnhsE  
           -0.0860395             -0.0819172             -0.0496835  
         CentralAir_Y        Functional_Maj2        Functional_Min1  
            0.0732900             -0.1285886              0.0054795  
      Functional_Min2         Functional_Mod         Functional_Sev  
            0.0247964              0.0169903             -0.3563441  
       Functional_Typ           PavedDrive_P           PavedDrive_Y  
            0.1103448              0.0005071              0.0463707  
SaleCondition_AdjLand   SaleCondition_Alloca   SaleCondition_Family  
            0.1887456              0.0061140             -0.0288674  
 SaleCondition_Normal  
            0.0385597  
# Resumen del ajuste
tidy(reg_fit_bic)
# A tibble: 55 × 5
   term        estimate std.error statistic   p.value
   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
 1 (Intercept)  11.8      0.0533     222.   0        
 2 LotArea       0.0354   0.00609      5.82 7.32e-  9
 3 OverallCond   0.0492   0.00455     10.8  3.37e- 26
 4 BsmtQual      0.0673   0.00556     12.1  3.80e- 32
 5 HeatingQC     0.0184   0.00487      3.77 1.67e-  4
 6 `2ndFlrSF`   -0.0548   0.00545    -10.1  4.38e- 23
 7 GrLivArea     0.196    0.00684     28.7  1.11e-142
 8 FireplaceQu   0.0267   0.00469      5.70 1.45e-  8
 9 GarageCars    0.0433   0.00555      7.79 1.31e- 14
10 GarageQual    0.0157   0.00430      3.66 2.66e-  4
# … with 45 more rows

Índices performáticos del modelo

reg_fit_bic |> extract_fit_engine() |>  performance()
# Indices of model performance

AIC       |      AICc |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------------------
-1627.362 | -1622.812 | -1331.336 | 0.888 |     0.884 | 0.133 | 0.136

Nuestro \(R^2\) es igual a 0.888, por lo que el ratio de información explicada del modelo es del 89 %. Este modelo es superior al ratio de información explicada con el modelo multivariante anterior.

Diagnosis

check_model(reg_fit_bic |> extract_fit_engine())

Supuesto de linealidad

En el gráfico, la línea es bastante atendencial, se ve bastante monótona en todo su recorrido (excepto en sus extremos). Comprobémoslo analíticamente con un ANOVA entre residuos y predictora:

adjustment <- 
  reg_fit_bic |>  extract_fit_engine()
lm(adjustment$residuals ~ adjustment$fitted.values) |>  anova()
Analysis of Variance Table

Response: adjustment$residuals
                           Df Sum Sq  Mean Sq F value Pr(>F)
adjustment$fitted.values    1  0.000 0.000000       0      1
Residuals                1458 25.971 0.017813               
lm(adjustment$residuals ~ I(adjustment$fitted.values^2) + adjustment$fitted.values) |>  anova()
Analysis of Variance Table

Response: adjustment$residuals
                                Df  Sum Sq  Mean Sq F value   Pr(>F)
I(adjustment$fitted.values^2)    1  0.0001 0.000061  0.0034 0.953411
adjustment$fitted.values         1  0.1195 0.119476  6.7338 0.009555
Residuals                     1457 25.8512 0.017743                 
                                
I(adjustment$fitted.values^2)   
adjustment$fitted.values      **
Residuals                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como se puede observar, el p-valor de la prueba F de Fisher-Snedecor es igual a 1, por lo que, a un nivel de significación del 0.05, no podemos rechazar la hipótesis nula, esto es, la no existencia de una relación lineal entre las predictoras y sus residuos. Por tanto, ello implica que no parece existir tendencia lineal entre residuos y predictoras en nuestro modelo. Lo mismo ocurre en el segundo caso: tampoco parece existir tendencia cuadrática entre residuos y predictoras.

Supuesto de homocedasticidad

En el gráfico, la línea sigue una fase decreciente al principio y creciente al final, se ve poco monótona en todo su recorrido. Comprobémoslo analíticamente con un test de heterocedasticidad:

check_heteroscedasticity(adjustment)
Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

Según el test de heterocedasticidad, el supuesto de homocedasticidad no se cumple para nuestro modelo. Esto, muy a menudo, resulta poco realista, pues el supuesto implica que la variabilidad del error de la variable objetivo es la misma para cualquier nivel de nuestra predictora.

ggplot(
  tibble("Observaciones" = 1:length(adjustment$residuals),
         "Residuos" = adjustment$residuals),
  aes(x = Observaciones, y = Residuos)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

En el gráfico, la recta de regresión se aprecia constante, en torno a 0, y los residuos se agrupan en torno a una banda también relativamente constante.

Supuesto de colinearidad

En el gráfico original de la diagnosis, prácticamente todas las observaciones están en la franja verde. Ninguna está en la franja roja.

Supuesto de normalidad de los residuos

En el gráfico, los percentiles empíricos de los residuos se acoplan bastante bien a la diagonal, excepto en los extremos. Comprobémoslo analíticamente con un contraste de normalidad:

library(olsrr)
ols_test_normality(adjustment$residuals)
-----------------------------------------------
       Test             Statistic       pvalue  
-----------------------------------------------
Shapiro-Wilk              0.939          0.0000 
Kolmogorov-Smirnov        0.0616         0.0000 
Cramer-von Mises         377.7789        0.0000 
Anderson-Darling         10.0147         0.0000 
-----------------------------------------------

Como se puede observar, los p-valor de los distintos test son iguales a 0, por lo que, a un nivel de significación del 0.05, podemos rechazar la hipótesis de normalidad. Por tanto, no podemos asumir a normalidad en nuestro modelo. Se han hecho pruebas con cambios en la receta añadiendo distintas transformaciones a nuestra objetivo (step_YeoJohnson(), step_BoxCox(), step_sqrt()), pero tampoco ha sido posible.

Influencia de los residuos

En el gráfico original de la diagnosis, todas las observaciones parecen estar contenidas entre las líneas punteadas al 0.5.

Fase 5.5: Evaluación del modelo (Regresión multivariante con selección de modelos)

Resumen

glance(reg_fit_bic)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik    AIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>
1     0.888         0.884 0.136      207.       0    54   870. -1627.
# … with 4 more variables: BIC <dbl>, deviance <dbl>,
#   df.residual <int>, nobs <int>

Nuestro \(R^2\) es igual a 0.8884, por lo que el ratio de información explicada del modelo es del 89 %.

Evaluación y predicciones en test

# Hacemos el split. Lo hacemos del 0.6 ya que hay pocos datos en train. 
set.seed(05492)

split_house <- 
  initial_split(house_train, prop = 0.6)

# Predecimos en test
reg_fit_bic_test <- 
  reg_fit_bic |> last_fit(split = split_house)

# Evaluamos en test
reg_fit_bic_test |>  collect_metrics()
# A tibble: 2 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard       0.149 Preprocessor1_Model1
2 rsq     standard       0.856 Preprocessor1_Model1

Visualización de errores en test

set.seed(05492)

# Errores en test
pred_test <-
  reg_fit_bic_test |>
  collect_predictions() |>
  mutate(error = SalePrice - .pred)
pred_test
# A tibble: 584 × 6
   id               .pred  .row SalePrice .config                error
   <chr>            <dbl> <int>     <dbl> <chr>                  <dbl>
 1 train/test split  12.3     2      12.1 Preprocessor1_Model1 -0.152 
 2 train/test split  12.2     3      12.3 Preprocessor1_Model1  0.0711
 3 train/test split  12.5     7      12.6 Preprocessor1_Model1  0.134 
 4 train/test split  12.3     8      12.2 Preprocessor1_Model1 -0.0565
 5 train/test split  11.7     9      11.8 Preprocessor1_Model1  0.0988
 6 train/test split  11.7    11      11.8 Preprocessor1_Model1  0.0742
 7 train/test split  11.7    13      11.9 Preprocessor1_Model1  0.224 
 8 train/test split  11.9    15      12.0 Preprocessor1_Model1  0.106 
 9 train/test split  11.9    17      11.9 Preprocessor1_Model1 -0.0135
10 train/test split  11.7    20      11.8 Preprocessor1_Model1  0.146 
# … with 574 more rows

En esta tabla se puede observar cómo, en nuestro modelo, las estimaciones se ajustan en gran medida a los valores reales de la objetivo SalePrice. Visualizaremos ahora estos valores a través de un gráfico.

g1 <- pred_test |> 
  ggplot(mapping = aes(x = .pred, y = SalePrice)) +
  geom_point(color = "#56BCC2", alpha = 0.6, size = 4) +
  geom_abline(intercept = 0, slope = 1, color = "#EB9891", size = 1.2) +
  theme_minimal() + 
  labs(title = "Resultados de la regresión lineal multivariante con selección de modelos",
       subtitle = "Los valores predichos deberían estar cercanos a la diagonal",
       x = "Predicciones",
       y = "Valores reales")

g2 <- pred_test |> 
  select(.pred, SalePrice) |> 
  gather(Distribución, value) |> 
  ggplot(aes(x = value, color = Distribución, fill = Distribución)) + 
  geom_density(alpha = 0.6) + 
  theme_minimal() + 
  labs(title = "Distribución de las predicciones sobre los valores reales de SalePrice",
       x = "Distribuciones",
       y = "Frecuencia")

multiplot(g1, g2)

Como se puede apreciar en los gráficos, las predicciones son bastante buenas, bastante mejores que en los anteriores modelos. A pesar de que se han incumplido algunos de los supuestos de regresión, la recta se ajusta bastante bien a los datos.

Comparación final de todos los modelos de regresión

Compararemos ahora los tres modelos de regresión: el modelo univariante, el modelo multivariante saturado y el modelo multivariante con selección de modelos.

compare_performance(reg_fit_uni |> extract_fit_engine(),
                    reg_fit_multi |> extract_fit_engine(),
                    reg_fit_bic_test |> extract_fit_engine())
# Comparison of Model Performance Indices

Name | Model |   AIC (weights) |  AICc (weights) |   BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------------------------------------------
..1  |    lm | -1932.0 (>.999) | -1931.9 (>.999) | -1916.1 (>.999) | 0.484 |     0.483 | 0.125 | 0.125
..2  |    lm |  -970.6 (<.001) |  -967.4 (<.001) |  -798.7 (<.001) | 0.891 |     0.887 | 0.133 | 0.136
..3  |    lm | -1030.6 (<.001) | -1023.4 (<.001) |  -772.7 (<.001) | 0.902 |     0.896 | 0.126 | 0.130

Como se puede observar, el mejor modelo es la regresión multivariante con selección de modelos. Presenta un \(R^2\) del 0.902 y un \(RMSE\) bastante bajo.

Predicción para el mejor modelo: regresión multivariante con selección de modelos

Nos quedamos finalmente con la regresión multivariante con selección de modelos (nos salió un 0.14492 de score :D).

# Predicciones
pred_bic <-
  predict(reg_fit_bic, house_test)
summary(pred_bic)

# Visualización de las predicciones para cada `Id`
final_pred <- 
  data.frame(Id = c(1461:2919), pred_bic) |> 
  dplyr::rename(SalePrice = .pred) |> 
  mutate(exp(SalePrice))
head(final_pred)

# Exportamos nuestro dataframe para subirlo a Kaggle
write.csv(final_pred, file = 'house_entrega_bic.csv', row.names = FALSE)